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Abstract 

Complex dynamical systems are often modeled as networks, with nodes representing dynamical 
units which interact through the network's links. Gene regulatory networks, responsible for the 
production of proteins inside a cell, are an example of system that can be described as a network of 
interacting genes. The behavior of a complex dynamical system is characterized by cooperativity 
of its units at various scales, leading to emergent dynamics which is related to the system's 
function. Among the key problems concerning complex systems is the issue of stability of their 
functioning, in relation to different internal and external interaction parameters. 

In this Thesis we study two-dimensional chaotic maps coupled through non-directed networks 
with different topologies. We use a non-symplectic coupling which involves a time delay in the 
interaction among the maps. We test the stability of network topologies through investigation 
of their collective motion, done by analyzing the departures from chaotic nature of the isolated 
units. The study is done on two network scales: (a) full-size networks (a computer generated 
scalefree tree and a tree with addition of cliques); (b) tree's characteristic sub-graph 4-star, as a 
tree's typical dynamical motif which captures its topology in smallest possible number of nodes 
and is suitable for time-delayed interaction. We study the dynamical relationship between these 
two network structures, examining the emergence of cooperativity on a large scale (trees) as a 
consequence of mesoscale dynamical patterns exhibited by the 4-star. 

We find a variety of coherent dynamical effects on the networks, which include: regular motion 
(emergent periodicity), weakly chaotic behavior (different from the uncoupled case), and self- 
organized motion characterized by close to zero Lyapunov exponents and anomalous diffusion in 
the phase space. Dynamical regions given as the intervals of coupling strength with distinctive 
motion and stability patterns are identified, suggesting a mesoscale interpretation of collective 
tree's dynamics in terms of 4-star's behavior. The system shows dynamical clustering in form of 
the structured phase space organization of orbits for all coupling strengths. Furthermore, various 
manifestations of the non-symplectic coupling are explored, including quasi-periodic orbits and 
strange attractors with weakly positive Lyapunov exponents. In our extended 4-star system 
whose dynamical units are driving each other, for certain coupling strengths we find the evidence 
of strange nonchaotic attractors displaying quantitative features which are known to appear in 
non-periodically driven maps. 

We employ the same two-dimensional chaotic maps for studying the stability of a real directed 
gene regulatory network of bacterium Escherichia Coli, the data on which are empirically known. 
The main cooperative effects including stability, clustering and long-range correlations are still 
present in the network's emergent dynamics. However, with increase of coupling strength the 
motion destabilizes on a sub- network of specific genes, although still maintaining some coherent 
properties. For comparison, a two-dimensional Hill model of gene interaction is implemented 
on the same Escherichia Coli network. We find the system to exhibit stable attractors and the 
flexibility of response to external stimuli, along with the robustness to fluctuations of the envi- 
ronmental inputs. 



Keywords: complex dynamical systems, gene regulatory networks, scalefree topology, modular 
networks, coupled chaotic maps, time delay, strange nonchaotic attractors, Escherichia Coli 
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Povzetek 

Kompleksne dinamicne sisteme pogosto modeliramo v obliki omrezij, kjer vozli predstavljajo 
dinamicne enote sistema, ki so sklopljene prek mreznih povezav. Primer taksnega sistema je 
omrezje genske regulacije, ki uravnava tvorbo proteinov znotraj celice. Tu vozli sistema ustrezajo 
posameznim genom. Znacilno za kompleksne dinamicne sisteme je sodelavanje gradnikov na 
vecih ravneh kompleksnosti, kar se odraza v dinamiki, ki je dolocena z ustreznimi funkcijami (oz. 
operacijami) sistema. Eden od glavnih problemov na podrocju kompleksnih sistemov je vprasanje 
stabilnosti delovanja sistema ob raznih notranjih in zunanjih motnjah. 

V disertaciji obravnavamo dvodimenzionalne kaoticne preslikave, skopljene prek neusmerjenih 
omrezij z razlicnimi topologijami. V sistemu uporabljamo nesimplekticno sklopitev, ki vkljucuje 
casovni zamik v sklopitvi med preslikavami. Stabilnost razlicnih mreznih topologij ocenimo z opa- 
zovanjem kolektivnega gibanja, ki ga dolocimo z analizo odmika od kaoticnega vedenja posameznih 
enot. Sistem podrobno obravnavamo na dveh ravneh: (a) na ravni celotnega omrezja (z numericno 
generiranim samopodobnim drevesom, ter z drevesom ki mu dodamo clique skupine); (b) na ravni 
drevesne mezostrukture, kjer obravnavamo podgraf '4-zvezda', ki je tipicen dinamicni vzorec za 
topologijo drevesa z najmanjsim moznim stevilom vozlov in je primeren za interakcije s casovnim 
zamikom. Proucujemo dinamicne odnose med tema strukturama, zlasti pojav kooperativnosti na 
ravni celotnega drevesa, kot posledico dinamicnih vzorcev ki ih 4-zvezda kaze na vmesni mezo- 
ravni. 

V raziskavi najdemo vpliv vec koherentnih dinamicnih pojavov: regularno gibanje (emergentna 
periodicnost), sibko kaoticno obnasanje (drugacno kot za nesklopljen primer), samo-organizirani 
razvoj s komaj pozitivnimi Lyapnovimi eksponenti ter anomalno difuzijo v faznem prostoru. V 
nekaterih obmocjih sklopitve se pojavijo znacilni vzorci razvoja in stabilizacije, ki sugerirajo na 
moznost mezoskopske razlage kolektivnega razvoja drevesa preko razvoja na nivoju 4-zvezde. 
Sistem kaze dinamicno zdruzevanje v obliki strukturirane organizacije orbit v faznem prostoru 
za vse vrednosti sklopitve. Poleg tega smo raziskali vpliv nesimplekticne sklopitve, ki se pokaze 
zlasti s pojavom kvaziperiodicnih orbit in cudnih atraktorjev s komaj pozitivnimi Lyapnovimi 
eksponenti. V sistemu razsirjene 4-zvezde, kjer so dinamicne enote sklopljene med sabo, smo 
za dolocene vrednosti sklopitve nasli cudne nekaoticne atraktorje, ki nakazujejo na kvantitativne 
lastnosti ki so znane za neperiodicno pogojene preslikave. 

Enake dvodimenzionalne kaoticne preslikave smo uporabili tudi za studij stabilnosti realnega 
omrezja genske regulacije na primeru bakterije Escherichia Coli, za katere so znani empiricni 
podatki. Dinamika omrezja v tem primeru se vedno kaze glavne kooperativne vplive, kot so sta- 
bilnost, zdruzevanje in korelacije dolgega dosega. Z narascajoco sklopitvijo se razvoj destabilizira 
enem podomrezju posameznih genov, kljub temu pa se dolocene koherentne lastnosti ohranijo. 
Za primerjavo smo implementirali dvodimenzionalni Hillov model genske interakcije na enakem 
omrezju Escherichie Coli. Pri sistemu opazimo stabilne atraktorje in fleksibilen oddziv na zunanje 
drazljaje, kot tudi neobcutljivost na motnje vplivov iz okolja. 



Kljucne besede: kompleksni dinamicni sistemi, omrezja genske regulacije, samopodobna (scale- 
free) topologija, modularna omrezja, sklopljene kaoticne preslikave, casovni zamik, cudni nekaoticni 
atraktoriji, Escherichia Coli 



Chapter 1 

Introduction 



This Thesis is a study of cohective dynamics of two-dimensional units 
interacting through the hnks of various network topologies. In this 
Chapter we introduce the basics necessary for understanding this 
work, which include the concept of complex systems realized through 
networks and the idea of network-coupling between the dynamical 
maps attached to the network nodes. We describe the process of gene 
regulation as the motivation behind the work presented here. In the 
last Section we precisely define the problem to be studied and set the 
directions of investigation. 

Nature presents us with an immense variety of complex systems which are composed of many 
simple units that through mutual interactions can achieve various patterns of coherent collective 
behavior [761 I82j . Their physically most intriguing aspect is that their behavior can often be 
very complicated and structured, despite the simplicity of their elementary building units. The 
self-organization effects are typically the key mechanism behind the cooperativity of the system: 
due to mutual interactions, system's units "communicate" and often manage to organize their 
functioning, with each single unit playing a particular role. Complex systems are characterized 
by their emergent behavior, which refers to the fact that the system's collective properties emerge 
as a consequence of interactions among the units, and are in general not exhibited by the isolated 
units without interaction. The nature of emergent behavior, along with the roles played by the 
system's units and small clusters of units are related to the origin of the system's complexity and 
its function. 

Examples of complex systems include climate, stock markets, telecommunication systems, 
World Wide Web etc [82] . All of these systems are inherently composed of many individual units 
(that can be atmospheric air masses, financial agents or web-pages), which develop cooperative 
behavior due to mutual communication/interaction, and/or due to being subjected to some ex- 
ternal influence. Most of the exploration in this field is focused on life and biological sciences, 
where perhaps the largest variety of complex systems can be found in [2T|. Social insects like 
ant colonies, exhibit a remarkably complex community structure with division and organization 
of labor that includes optimization procedures in searching for food [76]. Neural systems, widely 
investigated over the last few decades, explain the neural activity of a living organism based on the 
patterns of interaction and self-organization among the neurons, that are themselves very simple 
|13j . Ecological systems related to various environments, include very structured predator-prey 
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relationships that are sensitive and adaptive to different environmental conditions [21j . Gene 
expression regulatory system within living cells allows for efficient production of the proteins re- 
quired for functioning and development of the cell, based on the chemical inputs coming from cell's 
environment that signal the current needs of the cell [5]. Fish schools are able to quickly adapt 
to changes in environmental conditions like predator attacks, without any hierarchical structure 
or leaders [76] . 

The emergent function of a complex system is generally a non-equilibrium process, which 
depends on the external inputs and is sensitive to various parameters which are typically show- 
ing some sort of fluctuations. The functioning of gene regulatory system is, for instance, very 
sensitive to cellular inputs and characterized by the ability to quickly respond to them. Stock 
market quickly adapts to constantly changing actions of agents, as well as to different external 
circumstances. Functioning of complex systems with very different origins often share certain 
features that are more related to the needs of the optimal operation of the system, rather than to 
its particular purpose [82]. The operation performed by a biological system can differ according 
to the environmental conditions that come in form of external inputs; gene regulatory system 
produces different concentrations of proteins depending on the current needs of the cell [5] . 

There are various levels of interplay between the behavioral patterns displayed by the single 
units and the emergent behavior of the corresponding system involving many units. The nature 
of a complex system is typically either examined locally (on the scale of single units and their 
behavior), or globally (on the scale of the whole system). However, a complex system can also 
be studied at a mesoscopic scale, i.e. on the level of small clusters of units [40], which form 
functional modules. Modularity of biological systems is recently attracting a lot of attention, 
specifically in the context of gene regulatory systems where many functional modules have been 
revealed [51 [77]. Moduls are composed of few units which develop a certain level of cooperativity 
through particular interaction patterns. They play given roles in the emergent operation of the 
system, which is more involved than the behavior of a single unit, but still far simpler than the 
global system's functioning. The mesoscale approach to the complex systems bridges the gap 
between the single-unit level and the global level. 

Structural and organizational diversity of the complex systems are investigated through a vari- 
ety of disciplines ranging from experimental methods to mathematical and theoretical approaches. 
Besides studying natural systems, one is often interested in designing artificial complex systems 
that can serve as models for real (e.g. biological) systems. Computational modeling has a major 
role in this regard, as modern numerical techniques allow large and efficient simulations. By 
constructing the appropriate interactions among given isolated units, one can design a complex 
system with a specific technological or biological purpose. Due to numerical limitations, the iso- 
lated unit model is typically oversimplified; this however allows for big simulations involving very 
large number of units with possibly very complicated interactions [13]. The models are usually 
constructed with many parameters in order to allow for different types of collective behaviors 
to occur in relation to parameter variations. Last few decades saw a rapid development in this 
direction, particularly emphasizing the interdisciplinary nature of complex systems investigation 
[821 [5]. 

Computer generated complex systems can be designed in various ways. Typically, they can be 
represented as a network where a single unit is attached to every node. The interactions between 
the units are then assumed to follow the network links [9l [19]. In the context of real systems, 
the links between units are constructed following the directions of interaction among them, which 
results in a network of interacting units. In the case of artificially designed complex system, one 
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considers a pre-specified (large) network and attaches a simple (e.g. dynamical) unit with known 
behavior to every node. The cooperativity of the interacting units in then examined in relation 
to the network's architecture, nature of isolated units and external influences. Alternatively, a 
complex system can be constructed through a process of self- assembly: units are added one after 
another and linked according to pre-given rules producing a stable structure [21^ [82] . The system's 
emergent behavior evolves with addition of new units, and is not dependent only on the nature 
of isolated units and their interactions, but also on the properties of self-assembly process. This 
"bottom- up" growth is offer behind complex biological systems that arise through self-assembly 
of simple bio-units like cells or neurons [21] . 

In this Thesis we will study a complex system designed by attaching dynamical units to the 
nodes of a pre-specified network. 

1.1 Physics of Networks 

Networks are a typically used paradigm of representing complex systems, with nodes picturing 
individual units and network links specifying the directions and intensities of interaction [9l 119]. 
In a directed network links go exclusively from one node to another, whereas in a non-directed 
network each link symbolizes mutual interactions between the two nodes. Links can be weighted 
in case interactions between different pairs of nodes have different intensities. During the system's 
time-evolution, links and nodes can appear and/or disappear. 

Neural networks are ensembles of neurons interacting through synaptic contacts that can be 
represented as complex networks [21]. Predator-prey networks (food webs) are networks with 
species as nodes, and predator relations as directed links [3]. Telecommunication systems can be 
viewed as a network of users communicating through phone/Internet connections [H]. Regulation 
of gene expression involves networks of interacting genes that produce desired amounts of the ap- 
propriate proteins in response to the cell's needs, by directed activatory or inhibitory connections 
among the genes [5]- Social networks consider persons as nodes and friendship connections as 
links [Si [39]. 

Complex system realized through complex network topologies are often adaptive to certain 
inputs - by adjusting their interaction structure and/or topology according to set of environmental 
variables, they are able to modify their emergent operation and adapt to the assigned tasks. Study 
of adaptive networks is a fastly expanding field [37] • Adaptivity is recognized to be among the 
most prominent features of natural complex systems, which is relevant for their optimal design. 

The data behind the complex networks studied in physics usually come from experimental 
results found in certain measurement: World Wide Web (number of links per web-site), gene 
regulatory networks (directions of gene interactions), social networks (friendship relations) etc. 
The network's architecture in reconstructed upon these empirical data, and studied using various 
statistical and topological methods. Ultimately, many features of natural complex systems have 
been explored through studying architecture and statistical properties of the underlying networks 
[U [29] . Alternatively, as already mentioned, networks can be computer generated by given algo- 
rithms that produce connection patterns among a specified number of nodes. 

Elements of Graph Theory. The topological properties of networks are from mathematical 
viewpoint studied through graph theory \20\ I27j . We distinguish between the connected and 
unconnected graphs - the former is made of only one component (there is a path between any 
two nodes), whereas the latter is not. As in the case of networks, we distinguish between directed 



4 



CHAPTER 1. INTRODUCTION 



and non-directed graphs. Graphs are represented as ensembles of N nodes indexed by the index 
i = 1, . . . ,N with each node having ki hnks (to other nodes) which defines its degree. They are 
characterized through their topological properties, which besides the size (number of nodes) N, 
also include: 

- average degree < k > - total number of links divided by the total number of nodes < k >= 

- degree distribution P{ki) - the distribution of node's degrees, i.e. the probability that a given 
node has a certain number of links. It reports about global graph structure and is used to 
distinguish between types of graphs/networks (see below) |29| . 

- diameter - the maximal shortest path between two nodes belonging to the graph. Gives a rough 
idea of how compact is a graph: chains have biggest, and cliques (all nodes mutually linked) 
smallest diameters |27j . 

- clustering coefficient - characterizes how connected are a node's neighbors among them, which 
is often relevant for social (friendship) networks [8]. 

- modularity - expresses the structure of topological scales and presence of moduls in a network 

[iniETlEH] 

In computational simulations, a network/graph of size is usually defined through its x 
adjacency matrix Cij that specifies whether there exists a connection from node [i] to node [j]. The 
adjacency matrix is symmetric in the case of non-directed networks, while this is not necessarily 
true in the case of directed networks. 

In the work presented here we will deal with connected graphs/networks only, that will vary 
in size and degree distribution which shall be specified accordingly. We will be concerned with in- 
vestigations on local, global and mesoscopic scale, using both directed and non-directed networks. 

Networks Relevant for our Study. Different types of networks/graphs are referring to 
their architecture and/or their construction algorithms, which is particularly important for the 
applications in physics of complex systems. In the context of computer generated networks which 
are relevant for our study of complex dynamical systems, the considered types include: 

- Lattices and Chains - simple arrays of nodes connected deterministically in a form of a regular 
grid (e.g. linear chain with first-neighbor links) |27| . 

- Erdos-Renyi random graph - constructed by connecting a given ensemble of nodes by linking 
any two nodes with uniform probability [20j . 

- Small world networks (Watts-Strogatz model) - created from linear chains with random re- 
wiring of links between pairs of pre-selected nodes with a prescribed small probability |117j . 

- Scalefree networks - characterized by a power-law degree distribution (the probability of node 
having a certain degree is proportional to a power of that degree: P{k) ~ k~^), and grown by 
preferential attachment of the new nodes to the ones on the existing structures depending on their 
current degrees [29l |30] . 

- Modular networks - characterized by the presence of specific sub-networks occurring far more 
frequently than expected (with respect to the random networks), hence exhibiting a modular 
mesoscopic structure [51 [771 I78| . Modules (in dynamical terms also called motifs) of different 
sorts can be found in natural modular networks, having particular roles in network functioning. 

Natural and technological networks are widely studied in the recent years [HI [19] . Much have 
been understood regarding the topology and architecture of the networks that underline many 
complex systems, specifically in the context of relationship between the emergent function of 
a complex network and its topology |29[ [3]. In particular, many networks exhibit a scalefree 
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structure, with a power-law degree distributions having the power exponent between 7 = 2 and 
7 = 4. Scalefree structure was found in an ample and diverse variety of networks - World Wide 
Web, transport connections, social networks, gene regulatory networks etc. [H]. Some of these 
networks are natural while some other are technological, i.e. artificially designed: nevertheless, 
power-laws with very similar power exponents are found. It appears that hierarchical organiza- 
tion of links within these networks increases their efficiency in performing certain task [26] ; in a 
transport network for instance, hierarchy of connections characterized by existence of few "hubs" 
(well connected central nodes) reduces the path between any two nodes, thus allowing a more 
efficient transport [19]. Evolutionary advantages of scalefree topology were also emphasized and 
its role in the flexibility and operational robustness of many biological networks was underlined 
|84| . In addition, these networks often exhibit a modular structure, with identifiable families 
of motifs with different network functions [15^ I46j. In particular, gene regulatory networks are 
known to posses an architecture which is a combination of scalefree and modular topology [5] (to 
be described in detail later). 

Dynamics on Networks. Of particular interest for modeling complex systems is the study 
of dynamics on networks: isolated units (nodes) are described as simple dynamical systems in- 
teracting through the network links. Their interaction leads to emergence of collective dynamics, 
which has a specific relationship with the network topology, interaction parameters and the nature 
of the isolated units. The emergent behavior of neural or gene regulatory networks is primarily 
related to the dynamical features of single neurons /genes, and the properties of their interac- 
tion [211 US] • The efficiency of transport on complex networks is having an intricate relationship 
with the network structure [llOj . with the peculiarity of scalefree structure playing a crucial role 
[1081 1109] . Specifically, the transport of information in World Wide Web and similar complex 
networks is of great interest for a variety of modern applications [1061 1105j . The advantages of 
scalefree topology was moreover conffimed in a general context of network dynamics [84 115610.20] . 

Furthermore, study of dynamics on networks can give insights into the topology and structure 
of the networks, by observing the collective patterns and its respective transients |1241 110], The 
phenomena of synchronization (simultaneous motion of separate dynamical units) is of particular 
interest as it represents the most simple dynamical self-organization pattern [9]. General aspects 
of synchronization were investigated in a large variety of systems, including complex networks of 
dynamical systems |97j . interacting chaotic systems [89], systems of coupled cells [93], interacting 
strange non-chaotic attractors [94| and randomly coupled chaotic dynamical networks [72]. In- 
teractions between both continuous-time and discrete-time systems can lead to synchronization 
which is very easy to detect; in particular, more chaotic systems are also able to synchronize while 
maintaining their chaotic behavior [18j . 

Dynamics on network is studied in diverse areas of physics related to complex dynamical 
systems, including heat conduction in chains of alternate masses |63] . and in dynamic models of 
relaxor ferroelectrics [88] . 

The interplay between the complexity of emergent behavior and the dynamical nature of the 
isolated units can generally be of two types: 

• collective dynamics is very structured or possibly chaotic due to the nonlinear nature of 
interactions among its elementary units, that are themselves exhibiting very simple (linear) 
dynamics [57] . 
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• isolated units behave chaotically/erratically, but due to specific type of interactions among 
tliem, tlie emergent collective dynamics exhibits regular /periodic patterns with a high degree 
of internal dynamical organization |72| . 

In the present work we will examine a complex systems of the second type, with very chaotic 
individual units that exhibit a regular collective behavior due to neighbor-neighbor interactions 
through the links of the underlying network. 

1.2 Coupled Maps on Networks 

Discrete-time dynamical systems (maps) are often used for modeling complex dynamics, as they 
can exhibit a vast range of dynamical behaviors [49]. Physically, they typically represent dis- 
cretizations (Poincare maps) of realistic continuous-time systems, e.g. oscillators \119\ 1104] . A 
map is usually defined as iterative (recursive) sequence of values: 

Xn+i = f{xn), n G N or simply x' = f{x), x S (1-1) 

with function / determining the map's nature. Very simple maps often exhibit extremely rich 
dynamics that can be easily studied due to their computational simplicity. A famous example of 
one-dimensional map is the logistic map: 

f{x) = cx{l-x) (1.2) 

which displays a variety of behaviors when the parameter c is varied from to 4. This includes 
periodic oscillations, period-doubling cascade, chaotic intermittency and strong chaos (ergodicity) 
|104j . Logistic map belongs to the class of measure-preserving maps, whose dynamics conserves 
the phase space measure (linear length in this case), and is often used as a simple prototype of a 
chaotic system (displaying strong chaos for c = 4). 

A well-known example of two-dimensional chaotic system is the Chirikov- Taylor map, usually 
referred to as standard map [23], defined as: 

( ^' \ ua T\ f O + I + esinO mod 27r \ 

V/' )=^^^^^> = [l + esm9 ) ^^-^^ 

It arises as a Poincare section of the kicked rotator (frictionless plane rotation with a periodic 
kicking), with variables 6 and I denoting the angle and the angular momentum of the rotator, 
respectively. The map gives the values after the kick in function of the values before the kick. It is 
also a measure- preserving (area-preserving) map, as it is derived from a Hamiltonian dynamical 
system [231 166] . Parameter e that denotes the kicking strength is also the chaotic parameter 
for this map, in function of which the system Eq. (jl.3p displays an extremely wide spectrum of 
dynamical behaviors. They include periodic and quasi-periodic orbits, families of periodic islands, 
weakly and strongly chaotic behavior which has ergodic and mixing properties (demonstrated only 
numerically). As the standard map is a measure-preserving map that exhibits chaos, it represents 
a prototype of widely studied Hamiltonian chaos |66 1ll22j . Standard map will serve as an example 
of isolated unit used for designing the complex system studied in this Thesis. 

Network of Coupled Chaotic Maps (CCM) is constructed by assigning a map to each node 
of a given network and allowing them to interact through the network links due to dynamical 
coupling between them, giving a simple way to model dynamics on networks. The time-evolution 
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of a map/node on a network of CCM is given by the map's original update described by the 
Eq. (jl.ip . and an additional coupling contribution from its network neighbors at each iteration. 
The network is assumed to have N nodes/maps labeled by i = 1, . . . , A^. A coupled evolution 
equation for a CCM is usually written as: 



where the parameter controls the network coupling strength, and the sum goes over all the 
neighbors of the node i (where the normalization is done over the number of [i]'s neighbors 
denoted Ni) [68]. While the function / defines the original map, the function g represents the 
coupling (interaction) on the network that can be modeled in various ways. 

Network of CCM represent a simple and computationally easy model of complex systems, with 
logistic map being very frequently used as the dynamical unit (function / in Eq. (jl.4p ). Systems of 
this sort are extensively studied for various examples of underlying network topologies, interaction 
forms g and network coupling strengths /x, as they capture the essence of complex systems by ex- 
hibiting various cooperative phenomena |48^I49|. The essence of CCM network emergent behavior 
lies in dynamically correlated motion of maps at different network nodes, whose type of correlation 
generally depends on the network coupling strength. This primarily refers to the various synchro- 
nization effects [971 EH EZ] , including partial synchronization |123j , chaotic synchronization [181 E] 
and examples of entrainment of oscillators (that can also be modeled in continuous-time versions) 
[57\ I58j . Generalized synchronization phenomena include dynamical clustering [72\ 145] . driven 
phase-synchronization [H] and various other examples of dynamical self-organization [431 \6T\ 1102] . 
CCM were investigated for the cases of local and global coupling including external driving [44] . 
for regular lattices and complex networks [491 134] , with homogeneous and inhomogeneous cou- 
pling [121] and synchronous and asynchronous updating [1]. Recently, special emphasis was put 
on CCM collective effects on scalefree topologies [1241 168 ] [HT ] [60l [59] . small- world networks [35] 
and modular networks [101 183] , in the context of the mentioned developments in understanding ar- 
chitectures of these networks, including synchronization fittness of specific motifs [116] . Networks 
of CCM were also successful in modeling complex phenomena like dynamical phase transitions 
[2] and their classification [16]. Furthermore, coupled maps on networks are being increasingly 
applied for modeling of biological systems like cell-cycle dynamics [65] and regulation of gene 
expression [lllj . In this context the collective dynamical effects are known to be relatively robust 
to interaction parameters and perturbations [651 1107] . 

In regard of biological and technological applications of networks of CCM, recent investigations 
involve coupling between the units/nodes with a time delay. That is to say, a node receives input 
from the neighboring nodes "seeing" them in their dynamical states one (or few) time-steps in the 
past. Time delay in communication is ubiquitous in all natural systems. As no information travels 
instantaneously, real complex dynamical systems on networks ought to be modeled including a 
time delay in the node interactions. Synchronization and similar collective phenomena are actively 
investigated for time-delayed systems of CCM [73] for a vast range of topologies, interactions and 
lengths of time delay [69]. Specifically, for the systems of coupled logistic maps mentioned above, 
time delays are known to enhance the synchronization properties [111 l64] . 

However, the case of coupling between two-dimensional maps is still poorly investigated. The 
works that examine 2D coupled maps typically study the statistical system's properties only, like 
the anomalous diffusion [611 [6] , without investigating various non-symplectic coupling forms and 
their dynamical manifestations. As two-dimensional dynamics can be far more complex than its 
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one-dimensional counterpart, one expects that systems of CCM involving two-dimensional maps 
may have a larger spectrum of possible behaviors. Two-dimensional dynamics furthermore allows 
modeling of a wider range of complex systems, given that two degrees of freedom attached at each 
node give more coupling options, in addition to a richer local dynamics. Needless to mention, 
many real complex system posses elementary units whose internal isolated behavior cannot be 
well described by only a single degree of freedom (e.g. gene dynamics [Sj IllSj ). 

In this work we will be dealing with a time-delayed system of CCM involving two-dimensional 
maps, effectively showing the richness of its collective dynamical spectrum. As described in detail 
later, we examine CCM on large complex networks in comparison to the same dynamics on small 
sub-graphs (motifs). 

1.3 Models of Gene Regulatory Networks 

The regulation of gene expression is a process of fundamental importance for the functioning and 
growth of biological cells. Gene regulatory network represents a complex system of genes inter- 
acting by activating or repressing each other's expression level that defines the production rate 
of proteins [62]. Collective function of gene regulatory system responds to the cell's needs and 
provides the appropriate proteins at all times. Gene regulatory networks involve different scales 
of functions, and often have modular structure with specific groups of genes preforming particular 
tasks [281 E]. 

Process of Gene Regulation. Genes are building blocks of DNA that are responsible for 
synthesizing mRNA molecules from RNA polymerase present in the cell, through the process 
called transcription. The rate of transcription is determined by presence/absence of transcription 
factors, activators and repressors that can bind to the promotor region of a gene, located upstream 
from each gene on DNA strand that carries them. As pictured in Fig. II. H gene's transcriptional 
activity has three states depending on presence/absence of activators/repressors. Absence of both 
puts the gene in basal state (also termed leaky transcription) characterized by a faint transcription 
rate (much smaller than in active state) . Presence of an activator in the absence of all repressors 
defines active state of the gene when the transcription is maximal. When a repressor is bounded 
to its site, the gene is in inactive state regardless of the activators, and the transcription is 
fully inhibited. After transcription, mRNA initiates the production of proteins by the RNA 
translation process. Thus created proteins have various functions in the cell, including operating 
as activators/repressors for transcription of other genes [5| [62 } [TTS] . 

The operation of gene regulatory network is hence given by production of proteins by genes, 
that activate/repress other genes, creating collective functioning of the network that provides the 
products needed by the cell [5]. Besides proteins themselves, transcription factors can also be 
various environmental inputs able to modify the global operation of the network according to the 
cell's current needs. They are referred to as external transcription factors (ETF) |5|. 

Gene regulatory networks are directed networks: genes are attached to the nodes, while the 
regulatory interactions are pictured as the directed links. Moreover, they have two possible types 
of directed connections: activatory and repressory. Recent studies revealed their architecture 
in many details; in particular, for the cases of bacterium Escherichia Coli [112^ llOOj and yeast 
Saccharomyces Cerevisiae [113j where all of the data on gene regulation has been mapped. The 
network structure often shows scalefree properties, with few very connected genes (hubs) having 
central role in global network operation. In addition, these networks typically exhibit a modular 
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Figure 1.1: Process of gene regulation (with permission from jll8j ) with its three stages of 
operation pictured schematically: basal state, active state and inactive state. 



structure [771 S^i which has been extensively studied for the case of Escherichia Coli (E.Coli) 
[1001 [28| [96] . Modularity in the gene regulatory networks was found to be related to its global 
(biological) functioning [281 EI] and evolutionary design [TH [51] . Various computation algorithms 
have been designed for estimating the number and distribution of modules in a given network 
|52[ 199]. and the software applications are already available, like MFinder program by U. Alon [3]. 
The central issue of modular networks regards the relationship between the structure/ frequency 
of different modules and their functional role in the the network operation |115| . In the context 
of gene regulatory networks where modularity plays a crucial role, the presence of functional 
motifs has been revealed, whose topologies and network locations are directly related to the 
tasks they perform |28 [ 1100] . Specifically, the Feed Forward Loop (FFL) motif found in different 
regulatory networks is well investigated and its function fully understood, both theoretically and 
experimentally [5l[70]. As it appears, the coherent FFL always serves as a sign-sensitive delay 
element in all the gene regulatory networks investigated so far [71] . 

Appearance of modular networks is not limited to gene regulation systems only: functional 
motifs have been revealed in other real networks as well, notably in the case of brain networks 
[103] ■ metabolitic networks [461 [95] and some technological networks [77]. Various motifs can 
play different roles in networks of different origin, that may not be necessarily related to their 
topological properties, but to the needs and functioning of the network in question [92]. This was 
demonstrated on the example of Bi-Fan motif, another recurrent gene regulatory sub-network 
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(motif) [l2] . Modularity seems to have the central role in functioning of many natural networks 
|53| , indicating that a mesoscale approach might lead to applications in engineering of the systems 
that perform specifically assigned operations. Also, recent works are considering the possibility 
of engineering genetic circuits involving few genes that would perform pre-defined tasks, with 
potentially vast range of applications [32l [80l l4T] . 

Mathematical Models of Gene Interaction. As already described, genes interact by ac- 
tivating and repressing each others expression levels. They define the production rates of mRNA, 
which translate into proteins that may act as a transcription factors regulating the behavior of 
other genes [62] , Activation/repression processes are characterized by their respective thresholds: 
if the concentration of a given transcription factor is above a certain value, activation/repression 
starts [5]. A gene can be activated/repressed by other specific genes with interaction thresholds 
depending on the gene - its expression level will be determined by the interplay of these effects. 
The precise protein production rate is however not determined only by the quantity of activa- 
tor /repressor present (as long as this quantity is above the threshold), but by other interaction 
parameters. 

The gene interaction and the function of gene regulation network is known experimentally 
for some organisms, including those mentioned above \112\ \W0\ I96|. However, the mathematical 
modeling of gene interaction is not simple, as the details of this process can vary in time and from 
gene to gene in the same network. Still, a variety of models are currently in use, picturing the 
gene interaction with different levels of simplification [471 [85] . Among them are the following: 

• Boolean systems describe the interactions in truth-tables, where a gene at each time-step 
can be either "on" (activated and producing proteins) or "off' (repressed and not producing 
proteins), depending on the states of the genes influencing it [5l]. While in terms of dynamics 
this model is strongly oversimplified, it still offers a computationally easy model that allows 
large-scale numerical investigations of gene regulatory networks. Robustness and fiexibility 
in a boolean dynamical model of the E.Coli's gene network was studied recently [98]. 

• Step-function ID coupled maps are a more elaborate model, picturing threshold interactions 
as step- functions as follows [24]: 



where x[i] is the state of z-th gene/node, Tjj is the interaction threshold of node i influenc- 
ing node j, Kij are the interaction strengths and the binary value Sij (being either -1 or 
+1) describes the type of interaction (activatory/repressory). The step-function nature of 
activation/repression is captured in Heaviside function defined by: 



This is a simple and elegant model allowing more detailed studies of gene regulation, in- 
cluding theoretical studies of dynamical complexity and stability [671 I25j . In particular, 
robustness of gene regulation was recently confirmed in the context of time-delayed interac- 
tion using this model [55]. This model generalized the boolean approach, but also involves a 
simple summation of neighbor inputs on each gene similarly to the usual CCM model men- 
tioned above. A further generalization would include a truth-table for each gene, describing 
how neighbor inputs are processed for each gene separately. 
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Piecewise-linear ID Ordinary Differential Equations (DDEs) are the simplest continuous- 
time model of gene interaction constructed from the above discrete model |47] : 

= ^9ji{x[j]) - -nx[i], 
j 

where x[i] = x[i]{t) denotes the continuous-time state of the gene [i] at the time t. The 
piecewise-linear function gji describes the way gene [j] influences the gene [i], defining the 
input to the gene [i] by summing all the neighbor contributions, equivalently to the discrete 
model Eq. (jl.5p (again, instead of a sum, one could use truth-tables to specify relevance of 
different inputs for a given gene). The degradation constants ji describe the rate of self- 
degradation of each protein/transcription factor corresponding to gene [i]. This approach 
has also been amply studied both analytically [22] and computationally [311 150j . confirming 
a possibility of chaotic dynamics in gene interactions [[75] . It is also simple to implement 
numerically, although more directed towards concrete studies of dynamical phenomena, 
rather than collective network behavior. 

Hill-function 2D DDEs include both mRNA concentration (g) and protein concentration 
{p) as system variables, rendering internal gene dynamics two-dimensional, with protein 
production being governed by the concentration of mRNA. The equations are given by 
[SllIIH]: 

(1.6) 

dp[i\ p 

= "i^W - iiPiA 

where the function F^ (similarly to the function g above) models the gene to gene interac- 
tion, which is in this case either activation (+) or repression (— ), given by Hill's functions: 



(p) = f H ^ , activation 

F (p) = ^ + pn _^ j^n ' repression 



:i.7) 



This model has also been extensively studied theoretically [79j . and in particular, the two- 
gene system was analytically examined in detail for all combinations of interactions [118| . 
Note that this model besides being two-dimensional, also incorporates the largest variety 
of parameters describing the interaction: the Hill coefficient n determines the steep of ac- 
tivation/repression curve when the threshold T is crossed (at the limit of n — > oo function 
F^ becomes a step- function). We distinguish between non-cooperative binding (n = 1) 
and cooperative binding (n > 2) among the genes [118j . The factor /3 defines the maximal 
expression level that a gene interaction can produce. Both activation and repression expres- 
sions include leaky transcription state due to the presence of ^ > 0, which makes sure both 
F~^ and F~ are always bigger than zero. Note again that instead of summing over neighbor 
inputs, one can design more elaborate versions of input organization, as we shall describe 
in what follows. 
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This 2D continuous-time approach to the internal dynamics of a gene (termed Hill model from 
now on) will be used in this Thesis to design a model for gene regulatory network of E.Coli, which 
will be investigated along the lines of 2D network of coupled chaotic standard maps. It represents 
a two-dimensional model for gene interactions that locally exhibit regular dynamics, as opposed 
to the coupled maps system that locally displays only chaotic behavior. We will be confronting 
these two dynamical models through the same directed topology of gene regulatory network of 
E.Coli. The schematics of interactions in Hill model are shown in Fig. ll.2i 




Figure 1.2: The coupling structure of continuous-time 2D model of gene interactions [5l I118j 
showing local {q, p) variables at each node. They interact mutually and are influenced by external 
transcription factors (ETFs), with directed interactions being either activatory (red) or repressory 
(blue). 



Logic Gates. The remaining detail needed to complete the continuous-time model defined 
above regards the case when a gene is regulated by more than one transcription factor (among 
which can be p- variable of that gene itself, i.e. self-loop). The situation of more transcription 
factors regulating one gene is frequent, and the corresponding empirical data is usually provided 
gene by gene. As already mentioned, the simplest version of modeling this situation involves 
summation over all the neighbor inputs, while more detailed models include truth-tables or logic 
gates. Namely, some genes may need simultaneous presence of more activators to be maximally 
expressed, while some other genes enforce competition among many activators, and the final 
expression level is determined by the activator present in the highest concentration. The former 
case can be viewed as logic AND-gate, while the latter corresponds to the case of logic OR-gate. 
In naturally occurring regulatory networks logic gates in general have very complex operation 
patterns, and never operate as simple AND-gates or OR-gates. However, for the purposes of 
mathematical modeling we shall simplify the genes' logic gates, and determine the final expression 
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level of some gene regulated by N transcription factors pi through Hill functions Ff^ in one of the 
following ways: 

. OR-gate: QO = ma^j=i„„,N F^^{p[j]) 
• AND-gate: = mmj=i^,„^N F^^{p[j]) 

. SUM-gate: Gf = j^E^Zi F^Mj]) 

Although with lesser biological motivation, the SUM-gate has certain relevance, as it represents an 
analogy with the usual input-sum approach behind the typical design of coupled maps systems. In 
Chapter [5j we will be considering a dynamical model of E.Coli network, comparing two versions 
of the system with SUM-gates and AND-gates, as the logic gates involved in the E.Coli gene 
regulatory network appear to have properties of these two logic gate types [5]. 

1.4 The Subject of This Thesis 

The central subject of this Thesis is the computational investigation of collective dynamics of 
chaotic two-dimensional maps coupled with time delay on non-directed complex networks and 
their dynamical motifs. The Thesis summarizes the research results obtained in these directions 
and describes the computational techniques employed. 

Motivation. We study a network of coupled chaotic dynamical units, seeking the appearance 
of regular behavior due to the network interactions. Our aim is to test the ability of complex 
networks to induce stability into the dynamics of coupled chaotic units. In that context we consider 
two-dimensional standard map given by Eq. (ll.3p which exhibits strong chaos when isolated, and 
is characterized by phase space mixing and chaotic diffusion [66]. We develop a suitable model 
for testing the network stability by attaching a standard map to every node of a scalefree tree, 
which is taken as a representative of complex networks. In opposition to classical models involving 
one-dimensional units, we introduce two degrees of freedom at each node; many natural complex 
systems include units that cannot be appropriately described with a single degree of freedom. 
Gene interactions underlying gene regulatory networks are one of the examples, as they involve 
concentrations of proteins and mRNA [Sj. Our system also includes time delay in modeling 
communication among the dynamical units. 

We seek to recognize and study the departures from chaotic behavior of the isolated stan- 
dard maps generated by the network-coupling. Our hypothesis is that two-dimensional network- 
interacting systems generally posses equal propensity for self-organization as their one-dimensional 
counterparts, whose tendency towards cooperative behavior is well known [43 1 164^ [68]. 

As two-dimensional maps are far more dynamically rich, specially in the context of chaotic 
dynamics, our system provides insight into general nature of two-dimensional discrete-time oscil- 
lators with non-symplectic coupling. We study the behavior of separate nodes attached to the 
tree, and precisely examine departures from their nature as isolated nonlinear maps, as a further 
manifestation of the system's collectivity. 

Furthermore, we address the dynamical manifestations of network's modularity, by developing 
a mesoscale approach to understanding the emergent behavior of the system. We emphasize the 
role of the dynamical motifs, by examining the motion of the large scalefree tree in relation to 
the motion of its smallest typical sub- graph termed 4-star. The 4-star motif captures the tree 
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topology in only four nodes and is suitable for time-delayed interaction. We expect that the 
topological relationship between these two structures might be reflected into a relationship of 
their dynamical properties. This envisages the bottom-up approach to complex networks, where 
one seeks to reveal the properties of the large system by examining the properties of its local and 
mesoscale components - functional motifs. 

For comparison, we examine the standard maps coupled through the directed gene regulatory 
network of bacterium E.Coli. We claim that even a directed topology possesses equal tendency 
towards the dynamical regularity of coupled chaotic units. By investigating the stability of the 
emergent motion of standard maps on a real biological network, we are establishing a result that 
regards the general stability of the systems realized through the examined network. The result 
applies to other (less chaotic) nonlinear maps, demonstrating functional stability of the studied 
real network. As an extension, we examine the E.Coli's gene regulatory network model with 
two-dimensional (regular) units, showing its robustness and flexibility, which is expected by the 
biological origin of this network. We are showing E.Coli's topology to be able to generate self- 
organization that reflects the nature of the particular model, with both chaotic maps and regular 
Hill gene interaction model. 

The Statement of the Problem. We construct the system of Coupled Chaotic Maps 
(CCM) by considering a scalefree network (tree) and attaching a two-dimensional standard map 
given by Eq. (jl.3p to every network's node. The maps are allowed to interact through time-delayed 
coupling in the angle coordinate that goes along the network's (tree's) links. We comparatively 
examine the collective dynamics on the scalefree tree and a 4-node graph connected in a form of 
a star, representing the smallest dynamical structure capturing the tree topology. We study the 
emergent dynamical properties of this systems in relation to the variations of coupling strengths 
regulating the intensity of interactions among the network nodes/maps. We furthermore examine 
the dynamics of the same CCM realized through the directed graph of the gene regulatory network 
of E.Coli, and a two-dimensional Hill model of E.Coli's gene interactions given by Eqs. (jl.6p &: (jl.7p . 

Various dynamical regions of our systems of CCM are explored in relation to the coupling 
strength. The statistical properties of the emergent dynamics appearing on all networks are 
investigated using various methods, describing and quantifying the presence of collective effects in 
the systems. The stability of motion is addressed, demonstrating the ability of network structures 
to inhibit the chaoticity of isolated node's dynamics due to inter-node interactions. The dynamical 
relationship between the motion exhibited by the scalefree tree and its dynamical 4-star motif is 
investigated and discussed. We study the orbits of individual nodes attached to the network and 
investigate different phase space manifestations of non-symplectic interaction among the maps, 
like strange attractors and quasi-periodic orbits. A similar analysis is performed for the same 
CCM realized using the directed topology of E.Coli (which is downloaded from a database - see 
later). Relationships between the nature of dynamics on directed and non-directed topologies 
are discussed. In regard of two-dimensional gene regulation Hill model, we study the robustness 
and flexibility of the emergent behavior in relation to the environmental inputs and fluctuations 
(noise). 

The directions of our study and exposition of results will be the following: 

• appearance of regular behavior in network dynamics of coupled chaotic standard maps and 
the types of single-node orbits in relation to the coupling strength; dynamical regions with 
different types of collective dynamics for various ranges of coupling strengths for different 
network topologies 
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• general characterization of emergent dynamics, statistical properties of global motion de- 
pending on the dynamical region; transient times related to occurrence of regular behavior 

• stability of emergent dynamics and its relationship with dynamical regions; effects of non- 
symplectic coupling (strange attractors) and their characterization using non-linear dynam- 
ical systems approach 

• dynamical relationship between large structures (scalefree tree) and small structures (4-star 
motif) in function of the coupling strength; relationship between collective effects inducing 
emergent properties on large-scale vs. small-scale structure 

• dependence of emergent properties on the time delay, coupling form and standard map's 
chaotic parameter e 

• examination of the same CCM on a directed network of E.Coli: regularization properties, 
single-node orbits, dynamical regions, statistical characterization of the emergent motion 

• investigation of E.Coli's gene regulatory network designed via Hill model; system's robust- 
ness and flexibility of response to environmental inputs and noise 

Details and Organization of this Thesis. This Thesis is a summary of research results 
obtained in the directions mentioned above, performed as a requirement for Candidate's PhD 
degree. The research methods are purely computational, involving various programming codes 
written in C++ programming language, that compute the quantities of interest according to the 
investigation directions listed above. The obtained data are analyzed using smaller C++ codes, 
along with other software packages like MatLab and Python. 

As for the pictures representing the results, color plots were done with MatLab, while all other 
pictures were done using Gnuplot. The graphical representations of networks were produced 
in Pajek |17j . The research work was performed at the Department of Theoretical Physics at 
the Jozef Stefan Institute over a period of two and half years. The numerical simulations and 
data analysis were done using the Department's computing resources, specifically minos . ijs . si, 
eurus . i j s . si and zephyrus . i j s . si. 

The exposition of results in this Thesis follows the named directions, in form of Chapters 
dedicated to specific subjects. The structure of programming algorithms and the discussion of 
results is given along with the presentation. The Thesis is organized as follows: in Chapter [2j 
we introduce the structure of our network of CCM, along with the origin of different topologies 
to be examined. In Chapter [3l the general properties of CCM system are investigated in func- 
tion of the coupling strength. This includes main types of emergent orbits, dynamical regions 
and clustering, and the process of dynamical regularization. We also examine various statistical 
properties of the collective motion using network-averaged orbit approach. This Chapter finishes 
with a brief consideration of CCM system realized with other coupling forms. In Chapter HI we 
analyze in detail the stability of the emergent dynamics in relation to different dynamical regions, 
using three techniques: Finite-time Maximal Lyapunov Exponents, Standard Maximal Lyapunov 
Exponents and Parametric Instability Analysis. Occurrences of strange (nonchaotic) attractors 
and weak chaos, along with quasi-periodic orbits are found and examined. In Chapter O our 
CCM is investigated using the directed graph of the largest connected component of E.Coli gene 
regulatory network. A similar analysis is done in this case, discussing the differences with the 
case of non-directed topology. The same network is used for implementing a continuous-time 
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Hill model of gene regulatory network for E.Coli based on Eqs. (|1.6p fc (jl.7p . The expected stable 
behavior is found, and the flexibility of network's response to environmental changes and the 
robustness to noise is analyzed. The Chapter [6l is devoted to conclusions, including a systematic 
outline of the main results and the discussion of open questions. The extensive list of relevant 
references is provided in the end. 

Our Publications of Relevant for this Thesis. Most of the results presented in this 
Thesis have already been published. In particular: 

• Z. Levnajic and B. Tadic: "Dynamical Patterns in Scalefree Trees of Coupled 2D Chaotic 
Maps", ICCS 2007, Lecture Notes on Computer Science 4488, p.633-640, 2007 ([60]), con- 
tains preliminary results on scalefree tree collective dynamics from Chapter [3l 

• Z. Levnajic and B. Tadic: "Self-organization in Trees and Motifs of Two-Dimensional 
Chaotic Maps with Time Delay", Journal of Statistical Mechanics: Theory and Experi- 
ment, P03003, 2008 ([61]) contains the core material on statistical and stability properties 
of CCM on scalefree tree and 4-star reported in Chapters O andlH 

• Z. Levnajic: "Dynamical Regularization in Scalefree-trees of Coupled 2D Chaotic Maps", 
ICCS 2008, Lecture Notes on Computer Science 5102, p. 584-592, 2008 ([59]) presents in- 
vestigation of regularization of the collective dynamics on scalefree tree, covered in Chapter 

m 

• B. Tadic and Z. Levnajic: "Robust dynamical effects in traffic and chaotic maps on trees", 
Pramana Journal of Physics, 70, 6, p. 1099-1 108, 2008 ( |107j ) includes results on clustering 
of periodic orbits and dispersion of time series of CCM on scalefree tree, from Chapter [3l 

• Z. Levnajic and B. Tadic: "Collective Phenomena in Time-delayed Coupled Chaotic Maps 
on Directed E.Coli Gene Regulatory Network" is a work in preparation, aimed to include 
results from Chapter [5j regarding the CCM on directed networks 



The first two references from the list are provided at the end of the Thesis. 



Chapter 2 



Networks of Coupled Chaotic Maps 
with Time Delay 



The Network of Coupled Chaotic Maps (CCM) that we will study in this Thesis is designed 
by placing a Chirikov standard map Eq. (ll.3p on every node of a given network structure, and 
assuming the nodes/maps to interact via network links through a pre-defined coupling form. 

The Properties of Standard Map. For computation purposes we rewrite the standard 
map Eq. (ll.3p by introducing a new angular variable x = 6/271 and a new angular momentum 
variable y = I /2ii reducing the map's phase space to [0, 1] x M. The standard map now reads: 



with parameter e replacing the chaotic parameter e = 27re. We will use the standard map in this 
form to design our CCM, with prime (') denoting the map's update according to Eq. (j2.ip . 

As mentioned, the standard map represents a prototype of two-dimensional system displaying 
Hamiltonian chaos [66 ^ I122j . In particular, much attention has been devoted to its chaotic tran- 
sition as the paradigm of stochastic transition in dynamical systems \2>Q\ [38l I122j : at e = 0.971 
{e = 0.155) the phase space regions with locally chaotic behavior merge into a single zone of 
global chaos allowing transport through phase space. In Fig. l2.1l we show the phase space portrait 
of this map for two values of e displaying strong chaos. As standard map is an area-preserving 
system, remnants of the KAM surfaces persist at any e value, along with the portion of the 
phase space dominated by strongly chaotic dynamics (Fig. 12. Ik ). As in this work we are primarily 
interested in the effects that inter-node coupling has on regularization of chaotic dynamics, we 
shall focus on the standard map with the chaotic parameter fixed to e = 0.9 {e = 5.65). At this 
parameter e-value, the motion is strongly chaotic and mixing over the whole phase space without 



The chaotic properties of Chirikov standard map are discussed and 
the construction of our network of CCM is explained in detail. The 
coupling form follows the oscillatory nature of the standard map and 
models the interaction of real two-dimensional units. The network 
structures to be employed are exposed both graphically and through 
their topological properties. 




(2.1) 
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Figure 2.1: Phase space portraits for the standard map Eq. (j2.1|) done with 10 x 10 randomly 
chosen trajectories for e = 0.5 in (a) and for e = 0.9 in (b). 



any KAM surfaces remaining, as illustrated in Fig. 12. lb . Strong chaos at this parameter value is 
characterized by ESI [Ml M |38] : 

• normal diffusion in angular momentum coordinate y governed by < >= D^t with constant 
Dq for the regime of strong chaos given by Dq = e^/2 = 7r^e^/2 

• positive Kolmogorov-Sinai entropy hxs = = ln(7re) measuring the rate of informa- 
tion loss due to strong phase space mixing 

• positive Maximal Lyapunov Exponent (2D area-preserving system has at most one positive 
Lyapunov Exponent) which is equal to Kolmogorov-Sinai entropy by the Pesin Theorem 

The remarkable dynamical properties of standard map are a continuous source of past and recent 
explorations |101j . 

2.1 Construction of Our CCM on Networks 

We design our Coupled Map System by assigning a standard map Eq. ()2.ip with e = 0.9 to every 
node of a given network, which is assumed to have nodes indexed hj i = 1, . . . , N . Each node [i] 
becomes a 2D phase space, whose dynamical state at a discrete time t is denoted by (x[i]t, y[i]t). 
We construct the coupling by following the usual framework of diffusive phase-coupling in x- 
variable (angle) of the oscillatory standard map Eq. (12.ip . as done in [68l HSl SSI [73] . A one-step 
time delay in the difference of the angle coordinate values of the network neighbors is added. We 
also adopt the usual schematic approach to diffusive coupling: 



(1 — /i) X (isolated map-update) -|- /U x (coupling) 



(2.2) 
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as in [68] and other works. The fuh system is defined recursively, as fohows: 




) 



(2.3) 



in the sense that each node's update is a function of its own current state and the previous states 
of the neighboring nodes. Prime (') denotes the isolated map's next iterate as in Eq. (l2.ip . while 
t stays for the global discrete time. The parameter fi measures the network coupling strength, 
ki is node's [i] degree and /Cj denotes the network neighborhood of the node [i]. The update of 
each node is therefore the sum of a contribution given by its isolated standard map-update (the ' 
part) plus a coupling contribution given by the sum of differences between the node's phase- value 
and the phase-values of the neighboring nodes in the previous iteration, normalized by the node's 
degree. The sum of these two contributions is then implemented as pictured by Eq. (j2.2p . 

Note that time delay is realized here by comparing the updates of the isolated nodes (the con- 
sidered node [i] is an update ahead of its neighboring nodes in Eq. (|2.3p ). Some authors however 
define time delay by comparing the current iterate of the full system (time t) with its previous 
iterate (time t — 1) \73\ [TT\ I64j . For the purposes of our study, we will remain with the present 
definition (the CCM with alternative definitions of time delay were examined as well, without 
observing qualitative differences). 

Discussion of the Coupling Form. In designing our system, we followed the oscillatory 
origin of standard map to construct an appropriate coupling form that models the interactions 
among these 2D discrete oscillators. We couple the maps in the phase (action) variables only, in a 
way to make each node/oscillator "attempt" to adjust its phase to all of its neighboring oscillators' 
phases, but with receiving a time-delayed information about their dynamical phase states. The 
CCM system can therefore be seen as a system of one-dimensional phase oscillators (in action 
rr- variables), which are moreover internally coupled to a second variable (angular momentum y). 
Our system represents a generalization of well-studied examples of one-dimensional oscillators (in 
both discrete-time and continuous-time) on diverse types of lattices/networks [681 EH]- It also a 
generalization of typical diffusive systems (which are usually considered as chains and lattices of 
ID units [431 [2]). in the sense of employing two-dimensional units. 

As mentioned, time delay comes from the natural context of realistic physical interaction/ 
communication between the units of a complex system. Given that we are dealing with maps, 
time delay is easy to implement in the form of iteration difference. The imposed time delay is 
kept constant for all nodes (following [11] and in contrast to [73]). For simplicity, we consider only 
one time-step delay applied to the first network neighbors; however, there is no specific reason 
why the number of delay-iterations should match the network distance between the nodes ( |73j 
considers random time delay values). For comparison, we also show the results for the equivalent 
non-delayed system (see end of Chapter O). 

We are using the coupling scheme given by Eq. (j2.2p as this way the parameter fi measures 
the ratio between the isolated- map and the coupling contribution. Also, (1 — /i) factor in front 
of the ?/- variable inhibits the standard map's diffusion in the angular momentum [38], therefore 
allowing more interaction among the nodes. Moreover, we normalize the coupling contribution 
by the respective node's degree, as it is typically done [681 ESI E]) in order to make each node 
receive equal total amount of signal (and hence leaving the signal's two-dimensional direction to 
define the effect that the network neighborhood is having on the node at each time-step). 

This coupling form is non-symplectic, i.e., it does not maintain the area-preserving nature of 
the original isolated map, as opposed to the coupling form used in the context of coupled standard 
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maps in [6l[T2]. This comes from both time delay and (1 — /u) factor through which the coupUng 
is reahzed. Non-symplectic couphng is also typical in the context of CCM models of complex 
systems, actually it is rather difficult to construct a naturally motivated symplectic coupling, 
specially in the context of 2D coupled systems. Non-symplectic (dissipative) nature of coupling 
will actually enhance the synchronization and other collective proprieties of our system. 

Finally, as shown in [61] , the x-coordinate part of our system can be seen from the continuous- 
time point of view. By introducing a space-continuous variable x'{i,t) = x[i]'(t) (defined over the 
discrete network node-space and continuous time), the x-part of our CCM equation Eq. (|2.3p takes 
the form: 



x'{z,t + 1) = (1 - + f [E,eKS^[j]' - + Zje^S^j] " ^j]') 

(1 - + f \v'x[i]' - dtx[3\] - [(1 - /.) + ii(V2 - hdt)] x'{i,t) 



(2.4) 



which emphasizes the phase-diffusive nature of our network coupling (last expression in Eq. ([27 
and relates it to the usual expression for phase-coupling in CCM. 



2.2 The Network Structures Employed 

We study the collective dynamics of CCM Eq. (j2.3|) on networks with various topologies and N 
nodes. Computationally, we identify each network by its x adjacency matrix Cij. The 
results for different topologies are systematically compared. In the rest of this Section we expose 
topological details for each considered network and show their graphical representations. 



1. Scalefree Tree with N = 1000 nodes and symmetric links shown in Fig. l2.2[ This is the 
simplest example of the scalefree topology. The network is computer generated using the procedure 
of preferential attachment by 1 link/node [29], which is run until the size of = 1000 nodes is 
reached. We grow the network by adding a new node [i] at each step which is symmetrically 
attached to one of the existing nodes [j] with the probability 

L„(^n + a) 

P([j]) evolves with growth of the network and is re-calculated at each attachment-step. For the 
purposes of our study we fix a = 1. As it a tree structure, this network has no loops and there 
is a unique path between any pair of nodes. The scalefree tree's degree distribution is shown in 
Fig. 12. 7k in blue (averaged over 100 such trees), and has a clear power-law behavior: 

P{k) = const X fc-T. (2.6) 

The distribution's slope is 7 = 3 as expected, given the analytically known results 7 = 2 + a 
which holds for the case of scalefree trees [29] . In our investigations we always consider the same 
scalefree tree with N = 1000 nodes, grown as described above. 

The system of CCM is constructed from this network by assuming each node to represent a 
standard map which is coupled to its network neighbors in accordance with Eq. (j2.3p . Due to sym- 
metric (non-directed) tree's topology, every pair of connected nodes/maps influence each other. 
As mentioned above, our system of 2D CCM can be seen as a network of ID CCM with each 
map/node influencing itself (cf. Eq. (j2.ip ). which represents a dynamical self- loop (recall that tree 
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Figure 2.2: Non-directed computer generated scalefree tree with N = 1000 nodes. 



does not posses topological self-loops). Throughout this Thesis we will be primarily concerned 
with the emergent properties of CCM on this network. 

2. 4-star Motif which is a typical tree's sub-network with four nodes connected in a form of a 
star, shown in Fig. 12. 31 along with the schematics of coupling structure between the maps. This is 
the essential dynamical motif of the tree topology, capturing its structure in only four nodes and 
having symmetric (non-directed) coupling. It can also be thought of as the simplest non-trivial 
graph (smallest graph that is neither a chain or a clique). It displays graph diameter of 2 (links), 
necessary for introduction of one time-step delay into system of CCM (in a sense on time-delayed 
interplay between the central node and the branches). We will systematically examine 4-star 's 
dynamics in comparison to the tree's dynamics. We will also refer to the 4-node motif in the form 
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branch 




Figure 2.3: The 4-star motif with a schematic view of the coupUng interactions between the 
central node (hub) and a branch node, as defined by Eq. (|2.3p . 



of clique (all nodes mutually connected), called 4-clique, in order to emphasize the peculiarity of 
the tree topology with respect to other topologies. The dynamics of CCM on 4-star represents an 
eight-dimensional discrete-time nonlinear system. The manifestations of non-symplectic coupling 
will be explored in detail from the prospective of nonlinear dynamics in Chapter [H 

3. Modular Network generated from the scalefree tree by adding clique motifs, shown in 
Fig. 12. 41 We use a modular network whose modules are explicitly known in order to explore its 
relationship with tree and 4-star, as well as with its own modules (cliques). This network is grown 
by a preferential attachment of tree nodes combined with a non-preferential attachment of cliques, 
following the algorithm: 

• a certain number of tree nodes vtree is randomly selected with uniform probability from 
the interval [vmin,Vmax]- New nodes are then attached to the existing structure one after 
another. The attachments are done preferentially (as for the tree above): each new node [i] 
links to a selected node [j] on the existing structure with the probability -P([i]) = ^ {k"\-a) • 
The probabilities -P([j]) are updated with each attachment of a tree node. As before, we fix 
Q = 1 for all the networks 

• a clique of random size WcUque G [wmin-,Wmax] is selected and attached to the existing 
structure. The attachment is done non-preferentially: one node belonging to the clique is 
attached to a randomly selected node on the existing structure with a uniform probability. 
After the attachment of the clique, all /c^-s are updated accordingly with clique's size 

• the growth proceeds with db new selection of certa^in number of tree nodes Vtree- 

In all 

further attachments, both tree nodes and clique nodes are regarded as nodes belonging to 
the existing structure 

The network growth starts from the a tree with 10 nodes, after which the named two steps are 
repeated until the network size exceeds N = 1000 nodes. As previously, all the links are symmetric 
(non-directed). Note that in our modular network the moduls are integrated in the network, and 
not appended as on a " Christmas tree" . 
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For the purposes of our system of CCM we took Vmin = 5, Vmax = 10, Wmin = 4 and Wmax = 6, 
to obtain a network from Fig. 12.41 Since the average number of tree nodes attached is 7.5 and the 
average size of attached chques is 5, the ratio of tree nodes vs. chque nodes is 3:2. Our modular 
network on average has 600 tree nodes and 400 cUque nodes, implying it contains on average 80 
clique motifs. The asymptotic degree distribution for this modular structure (averaged over many 
equivalent realizations of this network) is shown in Fig. 12. 7k in red. It also exhibits a power-law 
degree distribution P{k) = const x k~'^ with somewhat steeper slope than in the case of the tree. 
As opposed to the scalefree tree exposed above, modular network contains topological loops, as 
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Figure 2.4: Non-directed computer generated modular network with modules indicated in black, 
constructed by repeatedly attaching a clique of size 4-6, combined with attachment of 5-10 tree 
nodes, together totaling = 1000 nodes. 



cliques include loops (in addition to dynamical self-loops on each node due to two-dimensional 
nature of dynamical units). The system of CCM was implemented on this network assigning a 
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standard map to each node and allowing them to interact according to Eq. (j2.3p . equally treating 
tree node maps and clique node maps. 

4. Gene Regulatory Network of Bacterium E.Coli, which was experimentally found in 2003 
[1001 [70] and whose adjacency matrix we downloaded from: 
http : //www . weizmann . ac . il/mcb/UriAlon/ 

We consider only the largest connected component of this directed network having N = 328 
nodes, shown in Fig. 12.51 This is our main applicational example, where we examine the collective 
dynamics on a real biological network. Both CCM Eq. (j2.3p and the 2D gene regulation Hill model 
Eqs. p.6p &: ()1.7p will be investigated using this network. For the system of CCM we will perform a 




Figure 2.5: The largest connected component of the directed gene regulatory network of bacterium 
E. Coli (interaction types, self-loops and names of the genes are not shown). 
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study similar to the study of non-directed network, in order to observe the differences introduced 
through a directed topology of a real network. For the 2D Hill model of gene regulation we will 
see how the same network can display properties of a biological complex system, like flexibility to 
environmental inputs. 

The full E.Coli gene regulatory network as available in the mentioned web-site, originally 
includes A'^ = 423 nodes which have their biological names as shown in Fig. 12. 61 For this network 
the rankings of nodes according to their in-degrees and out-degrees was computed and is shown 
in Fig. l2.7b . Some newer E.Coli databases include more genes, but are oriented towards biological 
applications rather than towards physical (network) studies. 
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Figure 2.6: Full E.Coli gene regulatory network including genes' biological names (without self- 
loops and the types of interactions), according to the data from |100t 170]. 



For the purposes of our investigations, we will focus only on the largest connected component 
of the E.Coli network from Fig. 12. 51 as it represent the dynamically most interesting part. Smaller 
connected components have their own dynamical behaviors which are unrelated to the dynamics 
of the largest component. In our exploration of the dynamics of CCM we will use the directed 
topology of E.Coli network only, while in the investigation of the 2D Hill gene regulation model 
the nature of gene-to-gene interactions (activatory or repressory) will also be included [1001 [70] 



26 



CHAPTER 2. NETWORKS OF COUPLED CHAOTIC MAPS WITH TIME DELAY 



(see Fig. 15. 171 in Chapter [5l for representation of the same network with types of regulatory 
interactions) . 




r 

(a) (b) 

Figure 2.7: Degree distributions averaged over many samples of the networks used for our systems 
of CCM with A'^ = 1000 nodes; scalefree tree like the one from Fig. 12. 21 (blue), and tree with 
addition of ~ 80 cliques such as the one in Fig. 12. 41 (red), fitted with respective slopes in (a). 
Ranking distribution of in-degrees and out-degrees for all nodes of full gene regulatory network of 
E.Coli's from Fig. 12. 61 fitted with q-exponential distribution (which is to be introduced in Eq. (l3.6p 
in (b). 



Chapter 3 



Collective Dynamical Effects in CCM 
on Networks 

Numerical implementation of our CCM on networks is described, and 
the main types of orbits exhibited by the system are shown. The 
regularization process leading to periodic orbits is explored. The 
dynamical regions related to various types of emergent dynamics are 
identified as intervals of coupling strength. A variety of collective 
effects arising due to inter-node interactions are recognized and 
statistically investigated, in particular the clustering of periodic orbits 
and return times with respect to phase space partitions. The Chapter 
ends with a brief discussion of networks of maps with alternative 
coupling forms. 

In this Chapter we identify and investigate examples of collective effects in dynamics of net- 
works of CCM defined in the previous Chapter. The study of our system is build upon observing 
its time-evolution for various values of coupling strength, using standard techniques of statistical 
characterization of motion, like return times distribution with respect to a phase space partition. 
We will focus on analyzing the departures from strongly chaotic properties of the isolated (un- 
coupled) standard map, which testify the presence of cooperativity among the system units. In 
particular, we will examine if the system's emergent behavior for small coupling strength can be 
characterized as regular, which would be a key argument towards the stability of real complex 
systems described through networks of this type. 

3.1 Numerical Implementation and Studied Types of Orbits 

We investigate our system computationally, by running numerical simulations in C++ programming 
language that implement the dynamics of Eq. (j2.3p with a pre-selected network topology. Each 
run consists in the following steps: 

1. set the network's size N, total number of iterations and the coupling strength /i 

2. introduce the topology by uploading the computer generated network adjacency matrix Cij, 
which has non-zero values Cij = 1 if a link between the nodes [i] and [j] exist 
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3. select the initial conditions for all the nodes by taking them randomly from (x, y) G [0,1] x 
[—1,1] with uniform probability 

4. iterate the network of CCM by updating all the x and y values for all the nodes simultane- 
ously, for the chosen total number of iterations 

5. analyze the final network state by extracting the information on the desired dynamical 
properties, which possibly includes application of additional program subroutines 

We examine the coupling strength interval of < /i < 0.08, focusing on certain regions more 
specifically as discussed below. As mentioned, standard map's chaotic parameter e is always fixed 
at e = 0.9 (regime of very strong chaos for the isolated standard map). The dynamics is typically 
run for 10^ iterations in order for transients to settle. The emergent dynamics describes the 
system's motion after the transients, which is also the final dynamical state of the system: once 
the transients are over, the system's overall motion pattern does not qualitatively evolve further. 

In our study, depending on the dynamical aspect under investigation, we examine system's 
orbits following three different approaches (to denotes the length of transients) , as in [61] : 

• orbit of an individual node [i] given simply as: 

(^[^]t,y[^]t)t>to' i = h---,N (3.1) 

which allows the usual study of nonlinear dynamics exhibited by a 2D discrete-time map (a 
node) coupled to other nodes, using standard dynamical system's tools. 

• network-averaged orbit (n.a.o.) defined as: 

1 ^ 

i^t,yt)t>tg = -^^{x[i\t,y[i\t)t>to- (3.2) 

i=l 

reduces the network's 2N-dimensional dynamics to the evolution of a single 2D phase space 
point. Useful in qualitatively illustrating global network's behavior, particularly in the 
context of statistics and stability of the dynamics 

• time-averaged orbit (t.a.o.), defined for each network node as: 

1 * 

(x[z],y[i]) = lim - — - V (x[i]fc, y[i]fc), (3.3) 

t^OO t — iQ 

k=to 

reduces the whole orbit of a node to a single phase space point that qualitatively captures 
the emergent motion. Good for examining the differences among the nodes in the network's 
final dynamical state. 

The employed programming codes are structured as follows: the main program starts by set- 
ting the simulation parameters according to steps 1., 2. and 3. from the algorithm above. A 
separate subroutine is called, which iterates the coupled dynamics for the desired number of iter- 
ations until final state is reached. A new subroutine is then called from the main program which 
performs the relevant analysis of the final dynamical state or simply prints one of the system's 
orbits as described above. 
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Orbits of a Single Node attached to Network. To illustrate the behavior of our network 
of CCM, we show some orbits exhibited by single nodes attached to a network (tree or 4-star), 
for various coupling strengths. In Fig. l3.1l we show typical (emergent) orbits exhibited by the 
branch-node (leaf-node) of the 4-star, which are also typical for many nodes on the scalefree tree 
for the same //-values. They are characteristic orbits corresponding to three main dynamical 
regions that we will examine in detail throughout this Thesis (the regions will be defined precisely 
as intervals of /x- values in the next Section) . The difference between dynamical regions is given by 
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Figure 3.1: Three typical emergent orbits appearing on the branch-node of the 4-star (a~c), with 
their y-coordinate time-evolutions (d-f). A chaotic orbit at ^ = 0.007 in (a&d), a periodic orbit 
at = 0.012 in (b&e) and a strange attractor at fi = 0.049 in (c&f). Inset in (c) zooms the left 
part of the orbit for a better visibility (we will often show only a half of an orbit). 



the general motion properties which are captured by these three orbits: chaotic motion persists 
at small /i-values (Figs. [3TTk fed), becoming regular/periodic with increase of /x (Figs. [3Tb & e), 
and exhibiting interesting dynamical phenomena for even larger /x, including strange attractors 
(Figs. [3TTb &:f). While standard map itself does not exhibit strange attractors (as it is an area- 
preserving system) jll9j . our network of CCM with non-symplectic (dissipative) coupling allows 
this possibility. 

Additional types of orbits are also displayed within the studied coupling strength range of 
< /i < 0.08 as shown in Fig. l3.2[ Quasi-periodic orbits Figs. l3.2b .fcd.b&:e construct another 
dynamical region, although the main attention shall be devoted to three types of motion mentioned 
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above. Strange attractors also appear on the 4-clique motif for the same couphng strength, as 
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Figure 3.2: Additional orbits of our CCM system. For simplicity, only a half of each orbit is 
shown (as in the inset of Fig. 13. lb ), together with its y{t) time-evolution. Quasi-periodic orbits 
on 4-star's branch node at /x = 0.014 (a&d) and at fi = 0.027 (b&e). Strange attractor on 
4-clique's node (all clique's nodes are topologically equivalent which is not the case with 4-star) 
at n = 0.049 (c&f). 



illustrated in Figs. 13. 2b & f. Clearly, our network of CCM posses a rich dynamics which is sensitive 
to both topology and the coupling strength. The (emergent) orbits on single-nodes of the scalefree 
tree are qualitatively similar to the orbits mentioned here (except for strange attractors), for the 
same /i-values. 

3.2 The Process of Dynamical Regularization 

In this Section we investigate the process of creation of regular orbits as a consequence of inter- 
node interactions during the time-evolution of the coupled dynamics, termed regularization pro- 
cess. We will quantitatively show how the chaotic dynamics of isolated standard maps becomes 
self-organized, by examining the increase of number of nodes exhibiting periodic orbits. 

Dynamical Localization and Periodic Orbits. The dynamics for very small coupling 
strengths is characterized by chaotic irregular behavior of all nodes, irrespectively of the underlying 
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topology (cf. Figs. l3.lb fc d). However, this chaotic motion does not exhibit normal diffusive 
behavior in angular momentum as in the case of isolated standard map [661138): instead, the motion 
remains localized in y-coordinate to a band of certain width, despite still being chaotic within 
the band. The network of CCM inhibits standard map's chaotic diffusion due to the interactions 
between the nodes, localizing the motion to a band in y-coordinate (Fig. 13. Ik illustrates what 
we intend as a band). The band's width furthermore shrinks with time-evolution, eventually 
becoming a periodic orbit after a transient which depends on the coupling strength value. 

We expose this quantitatively in Fig. 13. 3b where we examine the time-evolution of the mean- 
square distance < y'^{t) > for the 4-star's branch node, averaged over many initial conditions. 
Observe that even for very small coupling vales the diffusion is inhibited after a certain time, 
and a sub-diffusive behavior of < > emerges with slope 7 ~ 0. We moreover examine the 
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Figure 3.3: Dynamical localization of orbits for various coupling strengths: angular momentum 
diffusion < y^(t) > in (a) and average band width stretching rate < (ymax — ymin)(*) ^ ™ (^)- 
Plots are obtained by averaging over an ensemble of orbits of 4-star's branch node. 



time-evolution of the maximal band widths < (ymax — yniin)(^) measuring it for the same node 
averaged over many initial conditions, with results shown in Fig. 13. 3b . Again, for arbitrarily small 
coupling strengths the stretching rate is suppressed after a certain time (which is to be identified 
as the transient time). Once a band width reaches very small values (~ 0(1)), the orbit's band 
does not shrink further, but instead transforms into a regular (periodic) orbit as in Fig. 13. lb . As 
we shall clarify shortly, the chaotic orbits as the one in Fig. 13. lb are of transitory nature, as they 
gradually disappear with the system's time-evolution. As it appears from Fig. 13. 31 the transient 
times are independent from the initial conditions. 

We also examine the distribution of the band widths for the case of scalefree tree, considering 
the band widths of all tree's nodes. In Fig. 13. 41 we show the results for three values of coupling 
strength, averaged over initial conditions. The increase of coupling reduces the band widths 
towards smaller and more uniform values, eventually shrinking the motion of nodes to a narrow 
region in phase space with a band width ~ 0(1). 

As the 2D phase space domain available for nodes becomes smaller, the interaction among them 
becomes more efficient in terms of motion's collectivity. The accessible region for x-coordinate is 
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Figure 3.4: Distribution of band withs of all scalefree tree nodes averaged over many initial 
conditions, for three values of coupling strength. A transient of 10^ iterations was not included. 



always [0, 1], whereas a (1 — ^u) term is present in front of the momentum variable, which tends to 
reduce the accessible region in y-coordinate. Moreover, the dynamics in x-coordinate is coupled, 
with nodes receiving each other's inputs and adjusting to them. This interplay defines the system's 
time-evolution, which for small coupling strengths results in shrinking of the accessible band in 
y-coordinate. Finally, the shrinking of bands results in the creation of regularity in the coupled 
motion appearing in the form of periodic orbits of various periodicities, which are the central 
collective phenomena of our system of CCM Eq. (j2.3p . 

Periodic orbits are oscillatory and characterized by a period value; orbit shown in Fig. l3.1b 
has period value of four. The chaotic region is essentially a transient region, as the inter-node 
interaction eventually creates periodic orbits for any non-zero (but small) coupling. However, 
the transient times before the emergence of periodic orbits significantly increases for very small 
coupling strengths (for /i < 10~^ it is typically 0(10'') iterations). This regards both tree and 
4-star: after a transient that depends on the coupling strength, all the nodes exhibit periodic 
behavior similar to the one shown in Fig. 13. lb . The phase space locations and other characteriza- 
tion of periodic orbits will be extensively examined in the Sections to follow, along with the case 
of larger coupling strengths where the emergent orbits may not be periodic, but exhibit other 
self-organizational properties. 

Dynamical Regularization. For the case of 4-star all the nodes achieve periodic orbits 
simultaneously. However, in the case of tree with 1000 nodes, periodic motion firstly arises on 
certain nodes and then successively spreads with time-evolution, eventually reaching every network 
node. The increase of fraction of tree nodes having periodic orbits allows a direct examination of 
the regularization process. We illustrate this in Fig. 13. 5b where we show the increase of number of 
tree nodes exhibiting periodic orbits in function of time-steps, for few values of coupling strengths 
(following [59]). As observed already, for stronger inter-node interaction, the regularization of the 
whole tree proceeds much faster, and quickly reaches the final steady state with all the nodes 
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Figure 3.5: Regularization process for the scalefree tree. Number of periodic nodes in function 
of time for a single initial condition and various /i- values in (a), distribution of transient times tc 
averaged over many initial conditions, with log-normal fit for = 0.012. 



periodic. Interestingly, for some coupling strengths the process of dynamical regularization never 
(or extremely slowly) reaches the entire tree (cf. /i = 0.003 in Fig. l3.5h ). We also report the 
distribution (over many initial conditions) of transient times needed for the network to settle in 
the steady state for a fixed coupling strength of /i = 0.012 in Fig. 13. 5b . 

It is also instructive to consider the regularization process through the network structure: in 
Fig. 13. 61 we show four stages in tree's time-evolution by coloring periodic and non-periodic nodes 
differently at each stage. Outer less connected nodes seem to be somewhat faster in acquiring pe- 
riodic orbits; the regularization process also seems to follow the tree's structure by spreading from 
one regularized node to another. Nodes neighboring periodic nodes do seem to reach regularity 
somewhat more easily (as they appear to be inter-connected among them in all the pictures). 
However, other outer nodes seem to remain chaotic despite neighboring a regular node for a long 
time, and finally reach regularity together with some highly connected central nodes. 



3.3 Non-periodic Orbits and Identification of Dynamical Regions 

Periodic orbits are present in the collective dynamics of COM system Eq. (12.3p at all non-zero 
coupling values. The entire network either develops periodic orbits on all nodes, or on none of 
them: depending on the /i-value, there is a fraction of initial conditions yielding periodic orbits 
and the remaining fraction that never leads to periodic dynamics. This pattern is universal for 
all the studied network structures and it will be employed to define the dynamical regions for the 
emergent dynamics of our system of CCM. We will define the dynamical regions by considering 
what kind of emergent motion corresponds to those initial conditions that does not lead to periodic 
orbits. 

We define a periodic orbit of the node [i] by looking at its individual emergent orbit, as defined 
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Figure 3.6: Dynamical regularization process of the tree for ^ = 0.001: the periodic nodes (blue) 
and non-periodic nodes (yellow) at different stages of CCM on networks time-evolution for a single 



initial condition - after 2 x 10^, 4 x 10^, 10^ and 3 x lO'' iterations. 



^6 



in Eq. (j3.ip . for all the network nodes. We search for ti > to such that 



\x\i 



to 



< 6 and 



\xWto\ 



\y\iho\ 



(3.4) 



where denotes the end of transients and 6 accounts for numerical errors. In the simulations we 
generally use 6 = 10~^. If the value ti is found within some considered range of iterations the 
orbit is defined to be periodic, otherwise the orbit is assumed non-periodic (of course, the results 
may depend on (5- value and the range of iterations ti > Iq). The difference 



7T\t 



\i] =ti-to 



(3.5) 



is defined to be the period of the periodic orbit for the node [i]. 

The parameter we use to distinguish between the dynamical regions is the fraction of non- 
periodic orbits fr{np), that we examine in function of the coupling strength //. We obtain the 
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value of this fraction by averaging over initial conditions and the nodes of the network. The 
value of fr{np) reports how many nodes/initial conditions for a given /i- value have not developed 
periodic orbits after transients. Due to the finite length of transients, some initial conditions for 
the tree might only lead to periodic orbits on a part of the tree (cf. Fig. l3.6p . Fully periodic 
dynamics on network, given by fr{np) = can be understood as an analog of synchronization, 
occurring in other systems of CCM as described in the Chapter [TJ 

In Fig. 13. 71 we investigate the fraction of non-periodic orbits fr{np) for the [0,0.08] range of 
^-values, comparing the tree with the 4-star. After initial stabilization to fr{np) = occurring 
at ^ = 0.012, both systems destabilize again, but with the portions of non-periodic orbits that 
vary with ranges of /i- values: we identify these ranges as dynamical regions, according to the 

chaos quasi-periodic attractors 

periodic 



4-star — e- 
m ? SC tree q 
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Figure 3.7: The fraction of non-periodic orbits in function of the coupling strength /i for tree and 
4-star, averaged over 10 and 1000 initial conditions, respectively, with periods up to tt = 10000 
considered. Four main dynamical regions are roughly indicated through their respective peaks 
and ranges of //-values. We will refer to this plot as for definition of dynamical regions. 



scheme indicated in the picture. Note that these /x-regions of destabilization seem to be (at least 
qualitatively) common to all the investigated structures, which already suggests a dynamical 
relationship between them. In particular, we focus on four main dynamical regions, defined by 
their ranges of //-values and characterized by specific collective behaviors emerging for the case of 
initial conditions leading to non-periodic orbits. 

• Chaotic Orbits Region at small coupling values < /x < 0.012 is characterized by chaotic 
motion of all nodes (cf. Figs. l3.1b &: d) that eventually evolve into periodic orbits after a 
given transient time of diffusive behavior (the region is somewhat shorter for the case of 
4-star). This definition is given in accordance with the transient of 10^ iterations. 
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• Periodic Orbits Region between 0.012 ^ /U < 0.026 exhibits periodic orbits of all the nodes 
at all topologies (cf. Figs. 13. lb &: e). with no other types of emergent motion (after transients 
that are long enough). The details regarding structure and phase space properties of periodic 
orbits will be provided in the Sections to follow. 

• Quasi-periodic Orbits Region between 0.026 < /U < 0.04 is identified by quasi-periodic orbits 
(cf. Figs. 13. 2b fc e) emerging for some initial conditions for both structures (the remaining 
part of initial conditions leads to periodic orbits). Quasi-periodic orbits are displayed by all 
or none of the nodes, depending on the initial conditions, and represent final steady state 
of the coupled dynamics. They will be addressed in more detail in the Chapter HI 

• Strange Attractors Region for large coupling strengths 0.04 < /i < 0.065 is characterized by 
a variety of strange attractors appearing in the dynamics of 4-star for some fraction of initial 
conditions (cf. Figs. 13. lb &: f). For the tree, this region is characterized by localized noisy 
emergent orbits, that despite not being periodic, display a high degree of self-organization. 
This region will be studied in much more detail later. 

The dynamics with coupling strengths above /i = 0.065 is again characterized by periodic orbits 
for all initial conditions, and will be not investigated further here. Despite that precise /i-range 
of dynamical regions depends on the topology in question, the qualitative analogy between 4-star 
and tree is to be emphasized, suggesting the expected dynamical relationship between these two 
topologies that derives from their network structures. This relationship will be further examined 
in terms of statistical properties of emergent motion. 

The robustness of dynamical regions with respect to topology is further studied in Fig. 13. 8b 
where we compare the tree's profile from Fig. 13. 71 with the analogous profile for the modular scale- 
free network shown in Fig. 12.41 Note that the same dynamical regions can be equivalently defined 
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Figure 3.8: Fraction of non-periodic orbits in function of coupling strength for a single initial 
condition with periods up to vr = 10000 considered. Comparison of scalefree tree and modular 
scalefree network in (a); comparison of all structures for the transient of 2 x 10® iterations in (b). 



here as well, despite modular network being structured differently. Furthermore, in Fig. 13.8b we 
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compare the profile for all the mentioned networks, computed for a much longer transient of 
to = 2 X 10^ iterations. The profiles seem to maintain the same dynamical regions, except for the 
chaotic region which appears shorter due to more initial conditions reaching final periodic steady 
states due to longer transient. 

As opposed to classical studies of ID CCM mentioned earlier ([Ml [45l [43]) our network of 
CCM presents more intricate dynamical behaviors with respect to the various values of coupling 
strength. This is primarily given by the dynamical regions, which are not typical for ID CCM, and 
so far studied examples of 2D CCM. In the remainder of this Chapter we will statistically explore 
the nature of the emergent motion in order to gain better characterization of all dynamical regions. 
Since periodic orbits can be easily quantified, we will persist in differing among the dynamical 
regions by referring to the properties of the non-periodic motion in function of the coupling 
strength 

3.4 Dynamical Clustering of Emergent Orbits 

The nature of periodic orbits of system of CCM Eq. (j2.3|) is always given by oscillations between 
two groups of phase space locations, one on the left and the other on the right half of x-axis, 
just as in Fig. 13. lb . This is to say the period of each periodic orbit is (at least) a multiple of two 
(four, in the case of Fig. 13. lb ). However, periodic orbits of more distinct nodes can sometimes 
approximately share the same phase space locations (both left and right), without necessarily 
synchronizing or having the same period. This occurs for many initial conditions, /i-values and all 
topologies: some nodes group their final periodic motion by oscillating between same two phase 
space locations. All the nodes (orbits) that share the mentioned phase space locations for some 
initial conditions will be called a cluster of nodes (orbits). 

The global phase space organization of periodic orbits in our CCM on networks is characterized 
by each orbit/node belonging to one of the clusters (this holds for both periodic dynamical region 
and other dynamical regions, in case initial conditions lead to periodic orbits). To illustrate 
this graphically, we show ten iterations of all the tree nodes put on the same 2D phase space in 
Fig. l3.9h . for the case of tree with fi = 0.012. Oscillations are organized into 11 clusters, with each 
cluster being composed of two phase space locations, similarly to the orbit in Fig. 13. lb . Each of 
the 1000 tree nodes belongs to one of the clusters, as it oscillates between one location on the left 
and one on the right side of the phase space, separated by a distance in y-coordinate of about 1. 
While the x-coordinates of all the clusters' phase space locations nearly overlap, the y-coordinates 
are equally spaced on the y-axis. 

This dynamical clustering of nodes can be conveniently quantified by referring to each node's 
time-averaged orbit (which is basically a phase space point) defined in Eq. (13.3p and denoted by 
{x[i], y[i]). It is easy to recognize that time averages of orbits/nodes will be clustered, as the nodes 
sharing a cluster will have similar time-averaged orbits. Specifically, we shall examine the values 
of y[i], as through them the clusters can be easily defined, given the nearly uniform separation 
of clusters in y-coordinate. We report in Fig. 13. 9b the distribution of y[i]-values corresponding 
to the situation in Fig. l3.9h (but averaged over many initial conditions): well-defined peaks are 
in a clear analogy with Fig. l3.9h . suggesting a phase space symmetric organization of nodes into 
clusters with respect to y-axis. 

We also examine the distribution of y[z] -values for the whole range of coupling strengths, com- 
paring all the examined structures. The results are shown in Fig. 13.101 where for each /i- value we 
report the histogram of y[i]-values on the log-scale, thus making a 2D color histogram of values. 
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Figure 3.9: Final steady state of the tree for /i = 0.012 after a transient of 5 x 10^ iterations. Ten 
iterations on the same phase space for aU nodes in (a), and distribution of y[i]-values for all nodes 
computed over 1000 iterations and averaged over many initial conditions in (b). 



The difference between the profiles for a 4-star's node (Fig. 13. 10b ) and the profile for a clique's 
node (Fig. 13. 10b ) indicates outer less connected nodes have more tendency to organize into clus- 
ters, even for smaller values of coupling strengths. This is also visible by comparing the profiles of 
scalefree tree (Fig. l3.10b ) and modular network (Fig. l3.10l j) ~ tree's structure more rapidly devel- 
ops clustered organization of nodes at smaller couplings. The clear structural similarity of profiles 
for 4-star's branch node and tree derives from tree having many branch nodes (nodes having only 
one link to the rest of the network) , whereas modular network (despite its scalefree degree distri- 
bution) contains many cliques. The similarity between profiles of 4-star and tree, together with 
the similarity of profiles for 4-clique and modular network, clearly suggest a dynamics-topology 
relationship between a large network structure and its dynamical motifs. 

The clustering appears around the same coupling ^- value where the regularity transition occurs 
(/U = 0.012 with the transient of 10^ iterations), thus allowing a view of the regularization process 
from the clustering prospective. Furthermore, as it appears from the profiles in Fig. 13.101 the 
attractor and quasi-periodic dynamical regions seem to produce only faint departures from the 
very organized cluster profiles (visible for 4-star with fj, ~ 0.05). This is because all emergent 
motions (almost) entirely respect the observed dynamical clustering in terms of time-averaged 
orbits; nevertheless, their precise phase space locations may differ, as clear for the example of 
strange attractor in Fig. l3.1b and quasi-periodic orbits in Figs. l3.2b .fcb. 

We furthermore study the distribution of network-distances of the nodes sharing a cluster for 
a single final state on the tree, which is reported in Fig. 13. Ill together with distribution of tree's 
topological distances [60]. The prevailing distance of 2-links between the same-cluster nodes 
seems to be common for all the clusters, regardless of their respective y[i]-values. Interestingly, it 
is at smaller network-distances where the difference between same-cluster distances and topology- 
distances is the largest. The origin of this self-organization effect might lie in the time delay that 
induces local dynamical organization, that is however not present globally (as at larger network- 
distances the distribution follows the topological profile). Finally, in Fig. 13. 121 we examine the 
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Figure 3.10: 2D color histogram showing the y[i]-values in function of the couphng /i. 4-star's 
branch node in (a) and a chque's node in (b) averaged over many initial conditions. All the nodes 
for a single initial condition for tree in (c), and the modular network in (d). 



network locations of the clusters by showing the picture of the tree (cf. Fig. l2.2p with the same- 
cluster nodes sharing a color. The network distance of 2-links is indeed often visible for all of the 
five considered clusters. 

The dynamical clustering represents a clear collective effect present in the emergent dynamics 
of our CCM system, that was characterized using time-averaged orbit approach. Clustering of 
some sort is often found in ID CCM, although typically in the context of synchronization |45j . 
The clustering found here refers to a "weaker" form of synchronization involving only shared 
phase space locations, but nevertheless indicating a high degree of dynamical order \60\ I61j. 
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Figure 3.11: Distributions of the topological network-distances for the tree, and the network- 
distances between the nodes sharing a cluster, for four different clusters. A single final state of 
tree with /i = 0.017 was considered. 




Figure 3.12: The network locations of nodes sharing the same cluster (colored with the same 
color) for a single initial condition on the tree with = 0.02. 
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3.5 Statistical Properties of Dynamics of CCM on Networks 

The emergent behavior of CCM Eq. (j2.3p on networks is governed by collective effects arising as a 
consequence of inter-node interaction of isolated chaotic maps. In this Section we quantitatively 
characterize global collective effects using various statistical measures. 

In this context the node-averaged (or network-averaged) orbit defined by Eq. (j3.2p and denoted 
as {xt,yt) can be conveniently used, as it captures the global time-evolution of the whole CCM 
system in a single trajectory that evolves in a single phase space. To illustrate this we show 
in Fig. 13.131 the time-evolution of our CCM system through node-averaged orbit starting from a 
single initial condition. The formation of oscillatory behavior in a single cluster is visible, as ex- 
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Figure 3.13: 1000 iterations of the node-average orbit on the tree for a single initial condition at 
fi = 0.012. After iterations in (a), 5 x 10^ iterations in (b) and 5 x 10^ iterations in (c). 



pected from the known emergent cluster structure of the phase space at fi = 0.012 (cf. Fig. l3.9b ). 
The regular oscillatory behavior arises from an initial "cloud" after it splits and localizes, which 
leads to regularization, as captured here by the node-average orbit [59]. Due to large number of 
tree nodes, the final state of node-average orbit shown in Fig. l3.13t is very robust to the initial 
conditions. Because of the mentioned properties, we shall often refer to the node-average orbit in 
studying global aspects of our CCM on networks. 

Return Times with respect to Phase Space Partitioning. We examine the distributions 
of return times associated with phase space partitions (sub-divisions), defined as follows: consider 
a partition of the phase space (e.g. a rectangular grid), and a system's specific orbit (node- 
average orbit for instance) . Consider a single element of the phase space partition and observe the 
times At (number of iterations) between the consecutive visits that orbit makes to this partition 
element. The distribution of the time intervals At considered over all the partition elements for 
a time much longer than transients we call return times distribution with respect to a given phase 
space partition. Alternatively, return times distribution can by created by considering separate 
node's trajectories and averaging over them (although this is numerically far more demanding), 
or by simply considering the distribution for an orbits of a single-node attached to the tree. 
The investigation of return times distributions reveals long-range correlations of the motion in 
question, which report the presence of self-organization effects. 

In Fig. l3.14l we plot the return times distribution for node-averaged orbits for various topologies 
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and for two specific coupling strength values, representing periodic motion (/i = 0.012) and self- 
organized motion in the attractor dynamical region (/i = 0.051). In the periodic dynamical region 
(Fig. l3.1ib ) the short return times At < 100 appear most frequently, which is expected as the CCM 
system motion exhibits periodic orbits with relatively small period values for all the structures. 
Bigger return times in this region for the cases of small structures show constant profiles followed 
by exponential decays, whereas large networks show some correlations in form of power-laws. 
However, in the attractor region examined in Fig. 13.141 ) both tree and modular network display 
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Figure 3.14: Return times distributions of node-averaged orbits for various topologies, averaged 
over many initial conditions. Partition in form of rectangular grid with 1000 x 100000 elements 
covering {x,y) S [0, 1] x [—50, 50] was used. The case of /.i = 0.012 in (a) and fi = 0.051 in (b). 



power-laws indicating a presence of long-range correlations in the dynamics of network of CCM 
at this coupling strength. Moreover, node-average orbit of 4-star also shows a similar power-law 
behavior, in opposition to 4-clique which exhibits fully periodic behavior similar to what seen 
Fig. 13. 14^ . Mentioned power-law distribution is a sign of anomalous diffusion in phase space, 
which presents a sharp contradiction with the chaotic uncoupled case. The relationship between 
these profiles again suggests the collective dynamics of large scalefree tree to be generated on 
a smaller scale of a motif that captures its topology in only four nodes: the 4-star. Also, the 
profiles for tree and modular network are almost overlapping for both /i-values, which indicates 
4-star to represent modular network's collective dynamics to some extent as well. This could be 
understood by observing our modular network to be constructed (cf. Chapter El) with many free 
links which become branch nodes, similar to the ones of the tree/4-star. 

In Fig. l3.15l i we show in more details the return times distributions for various dynamical 
regions for node-average orbit of 4-star, obtained by averaging over many initial conditions |61j . 
Chaotic behavior for small coupling strengths at the limit produces a constant distribution, as 
all return times are equally possible (this profile would however change due to regularization 
process, in the case when transients longer than 10^ are considered). In periodic region only 
small return times are prominent, while larger ones show similar behavior as in Fig. 13. 14^ . In the 
attractor region, the profile exhibits a structure similar to Fig. 13. 141 3. but without a clear slope as 
it comprises many initial conditions leading to different phase space structures. We examine the 
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Figure 3.15: Return times distributions and attractors for 4-star. Node-average orbit in three 
dynamical regions, averaged over many initial conditions in (a), two single-node attractor orbits 
occurring for // = 0.048 at the 4-star 's branch node, whose parts are shown in (c) & (d). Same 
grid partition was used, covering only the portion of phase space containing attractors. 



attractor region at /i = 0.048 in Fig. lS.lST p where we show the return times distributions for two 
single-node orbits (attractors) shown in Figs. 13. 15b &: d. Both attractors are exhibited by 4-star's 
branch node at this /x-value, but for different initial conditions. The return times distributions 
can be fitted with the q-exponential function, defined by |114j 

eix) = B, (^1 - (1 - q)^y'\ (3.6) 

with q = 1.6 for the attractor in Fig. 13. 151 :; and q = 1.5 for the attractor in Fig. l3.15B . We shall 
study these attractors further in the context of dynamical stability in Chapter HI 

Finally, we can distinguish between various dynamical regions by examining slopes of return 
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times distributions. We consider for simplicity only the j;-coordinate of the node-averaged orbit of 
the tree. The distributions for the periodic and quasi-periodic regions averaged over many initial 
conditions are reported in Fig. 13.16^ : note the similarity between two profiles referring to the 
periodic region at /i = 0.021 and /i = 0.034, and between other two profiles referring to the quasi- 
periodic region at /x = 0.014 and /i = 0.027. We also show in Fig. 13. 16b the distributions that 
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Figure 3.16: Return times distributions of tree's node-averaged orbit with respect to x-coordinate 
partition involving 100000 elements, averaged over many initial conditions, and for various fi- 
values. Transients were not included. 



regard the case of initial regularization at /i = 0.001, = 0.004 and the attractor region for /x = 
0.051. Each region seems to have a unique power-law slope in the return times distribution, as the 
dynamics of the tree appears to be dependent on the /x- value, with various types of collective effects 
present in all examined //-values. Although tree has many non-periodic orbits for all examined 
//-values, the presence of anomalous diffusion in the phase space motion indicates departure from 
the uncoupled dynamical regime. 

The patterns of power-law and q-exponential fit observed in the return times distributions for 
various topologies and coupling strengths, testify clearly about the presence of network collectiv- 
ity in the emergent dynamics of our system of CCM. We noted the analogies between 4-star's and 
tree's emergent statistical properties, pointing again towards a presence of a dynamical relation- 
ship between them induced by their known topological relationship. Similar statistical approaches, 
e.g. recurrence time statistics, are extensively in use in the context of coupled standard maps [6]. 



Period Values Statistics. System's periodic orbits emerging in periodic and other dynamical 
regions, can be described by both their cluster-structure organization in the phase space, and by 
their period values. For the period value associated with a given periodic orbit of the node [i] 
we intend the value 7r[i] as defined by Eq. p.Sp . The period of an orbit clearly depends on the 
choice of 5- value in Eq. (j3.4p : throughout this Section we will maintain 6 = as previously. We 
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checked our main results using different values of 6 (10~^, 10~^, 10"^'^), concluding their general 
properties to be qualitatively independent from this choice. 

In Fig. l3.17h we show the histogram of period values obtained for a single final periodic state 
of the tree at /u = 0.012. In the picture we include the periods up to vr = 3000 (some nodes have 
bigger periods) showing the number of nodes (out of total 1000) having certain vr-value. The 
period values appear in a discrete sequence, with (almost) each node having vr-value multiple of 
240, as it can be seen from the profile. This structure of values is robust to the initial conditions 
for this value of coupling strength. The discretization of periods is also found for other //-values, 
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Figure 3.17: Period values of tree nodes at a final periodic state. Histogram of period values (up 
to TT = 3000) for a single initial condition at fi = 0.012 in (a), distribution of period values for all 
nodes averaged over many initial condition for two /U- values fitted with a slope of -2 in (b). 



but instead of 240 the basic period (basic multiplier) can be different (smaller for larger coupling 
values). The property is moreover independent of (5- value in the definition of periodicity Eq. (|3.4p . 
although the basic period for the same coupling strength may vary. In Fig. 13. 17b we report two 
distributions of period values obtained by averaging over many initial conditions, for = 0.012 
(periodic region) and fj, = 0.08 (also a periodic region, although not specifically studied here, cf. 
Fig. l3.7p . Both distributions in the limit of large vr- values exhibit power-laws with slopes close to 
-2, again indicating a presence of long-range correlations in the final dynamical states at those 
/i- values. As expected, the case = 0.08 displays predominantly smaller periods up to vr = 10^, 
whereas for /i = 0.012 the periods range up to 10*^ (which was the maximal vr-value we measured). 
This effect is systematically present in our system - smaller interaction strength leads to richer 
collective behavior. The statistics of period values is difficult to produce in the context of 4-star 
as the period values are much smaller. 

The node-average orbit approach is also less convenient in the case of period values statistics, 
as it typically has extremely large periodicity (at least equal to the largest node period). However, 
it is instructive to examine the ranking statistics of period values for specific nodes in relation to 
the coupling strength, which can be easily done for all the topologies. The ranking is obtained 
by ordering periods from the largest to the smallest observed value, and hence associating a rank 
r with each period value vr, without repeating the period values. In Fig. 13. 181 we consider the 
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rankings of period values of 4-star's branch node with an outer node on the tree (a node far 
away from the hub with degree 1), obtained for many initial conditions for two coupling values. 
Again, the network of CCM with = 0.08 (Fig. 13. 18b ) mainly exhibits smaller period values, 
while the profile is broader for /_f = 0.012 (Fig. 13. 18b ). It is also clear the profiles resemble each 
other for both coupling values, in terms of shape, structure and the density of ranks with respect 
to the period value. This is expected from the network-location analogy between 4-star's branch 
node and the outer tree's node, and provides a further confirmation of the dynamical relationship 
between these two structures. 

It is to be noted periodic orbits are not common in the context of ID CCM. They are a 
distinctive feature of two- (and more) dimensional CCM which indicate collectivity of system's 
motion. Periodic orbits are very rare for the uncoupled standard map, which proves them to be a 
consequence of (non-symplectic) coupling in our CCM [59]. As opposed to other studies of CCM 
on lattices and networks (cf. Chapter [TJ), our system does not exhibit synchronized behavior. 
Nevertheless, given our context of 2D coupled units, the presence of dynamical clustering and 
emergent periodicity are a sign of self-organization in the system. 

The Dispersion of Time-series. Another illustrative way to study the collective properties 
of network of CCM on the tree is to look at the dispersion of the system's time signals, i.e. to 
investigate the correlation between signal's time average and its standard deviation. For the the 
case of our system we consider time series of angular momentum y[i]t as the signal, and compare 
their time averages y[i] with their standard deviations a[i] taken over the averaging interval. The 
results are reported in Fig. 13. 191 where we examine y[i] time series dispersion on two scales for 
various /i-values. In Fig. l3.19b i we show a log-log scale plot for all the tree nodes: for each coupling 
strength the variations of a[i] are very small in relation to the variations of |, in contrast to the 
general scaling law referring to the dispersion of time signals in (see |107j and references therein) : 



a[i] oc const x (y[i]) 



(3.7) 
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and typically displays exponents > 2- This profile also offers an alternative view on regular 
vs. chaotic dynamical region as the periodic clusters have a[i] < 1, whereas chaotic nodes show 
a[i] ~ 10. To illustrate this further we show in Fig. 13.191 3 the zoom of the cluster region of disper- 
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Figure 3.19: Dispersion of time series y[i]t, obtained by computing y[i] and a[i] over 10000 itera- 
tions for each tree node at a single steady state for various coupling strengths. Log- log scale plot 
for all nodes in (a) , and a zoom of the central part of the dispersion profile on linear scale in (b) . 



sion plot on a linear scale, with cluster organization of y[i] clearly visible for all coupling strengths. 
Observe the correlation between node's y[i]-value and their standard deviation values cj[i], which 
goes from no organization (// = 0.001), via clusterization process (/x = 0.007 and fi = 0.012), to 
fully clustered phase space organization (^u = 0.051). 

In this Section we have shown various statistical feature our system's behavior, which uniformly 
suggest presence of cooperativity in the emergent motion. All the found properties are in a clear 
contradiction with the case of uncoupled standard map. Furthermore, we found more evidence 
of dynamical relationship between tree and 4-star, which was initially known to exist only in the 
sense of their topologies. 

3.6 Networks of Maps with Other CoupUng Forms 

The coupled discrete-time equations defining our system are given in Eq. (j2.3p . They include a 
specifically selected time-delayed non-symplectic coupling form, and involve a fixed the standard 
map parameter e = 0.9. As described previously, through this particular coupling form we seeked 
to model a realistic coupled system of 2D oscillator units that would resemble a system of in- 
teracting genes. The particular system's set-up is also used to test the stability of networks of 
CCM realized through various network architectures. However, it is clear that Eq. (12.3p does not 
exhaust the possible versions of coupling forms, in terms of their structure and parameters. In 
this Section we will examine the CCM established through the same networks, but with different 
coupling forms and interaction parameters, showing the peculiarity of our specific choice. 
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The CCM without Time Delay. Our network of CCM is designed including the time delay 
of information traveling between the network's nodes that is ubiquitous for all natural systems. 
For the sake of completeness, we shall here present some results concerning the version of our 
system without the time delay. The coupled non-delayed equations equivalent to Eq. (j2.3p read: 



(l-A^) 



X l\ 



+ 







x\i\ 



(3.8) 



and in contrast to Eq. ()2.3p have prime (') on both x terms in the coupling expression, implying 
no time delay (an alternative version of non-delayed equations includes no primes at any of x 
terms). We will examine the collective dynamics of system of CCM defined by Eq. (l3.8p in a way 
analogous to previously examined time-delayed CCM, using the same network topologies. 

In Fig. l3.20l we compare the fractions of non-periodic orbits for delayed vs. non-delayed case, by 
examining scalefree tree and 4-star. Non-delayed case exhibits similar de-stabilization of emergent 
dynamics, with apparent dynamical regions. Although the profiles are similar in terms of their 
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Figure 3.20: Fraction of non-periodic orbits in function of the coupling strength /x, comparing 
delayed vs. non-delayed versions of network of CCM. Averaged over many initial conditions for 
4-star in (a), and for a single initial condition for nodes of the tree in (b). 



general shapes, the delayed case still presents somewhat wider dynamical regions in the sense of 
the corresponding coupling strength intervals. We furthermore investigate the statistics of node's 
time averages y[i] in Fig. l3.211 for non-delayed CCM of the 4-star and tree. Note that clustering 
regime is again present, but the persistence and the constancy of clusters is weaker, specially for 
larger /i-values in the attractor dynamical region. This indicates the non-delayed case exhibits 
less self-organizational properties than its time-delayed version (cf. Fig l3.10p . particularly in the 
attractor region which is most interesting in terms of collective effects. 

The motion of the non-delayed CCM does achieve regularity and periodic orbits organized 
similarly into dynamical regions, but does not posses the dynamical diversity of the time-delayed 
case. This suggests further explorations of time-delayed systems is needed, as they might present 
more realistic collective dynamical effects, which is expected given the universality of the time 
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(a) (b) 

Figure 3.21: 2D color histogram showing the y[i]-values in function of the couphng ^ for non- 
delayed CCM Eq. (l3.8p . done equivalently to Fig. 13.101 The 4-star's branch node over many initial 
conditions in (a), and a single initial condition for all nodes of the tree in (b). 



delay in nature [llj . 

The CCM with Different e-values. Throughout our investigations we keep the standard 
map's chaotic parameter fixed to e = 0.9 in order to be able to comparatively study different 
topologies and coupling strengths. It is necessary to consider a large e-value as we are interested 
in the interaction of strongly chaotic maps, examining the ability of various topologies to create 
regular motion through the inter-node interactions (standard map's chaotic transition occurs for 
e > 1.55 [Ml [36]). 

Here we investigate the system of CCM given by Eq. (j2.3p . taking different e- values on the 
scalefree tree, by considering the profiles for fractions of non-periodic orbits as done previously. 
The results are shown in Fig. 13.221 where we compare the systems with e = 0.85,0.95 and 1.0 
with already examined case of e = 0.9. While the network of CCM with e = 0.85 stabilizes very 
quickly, the system with e = 0.95 stabilizes only for a very large coupling strength (/i = 0.07), 
and the CCM system with e = 1.0 never stabilizes in the studied coupling range (although, it is 
interesting to note the e = 1.0 case exhibits regularity in a part of the tree for /u > 0.02). The case 
of e = 0.85 as the only one that achieves regularity, shows very small portions of initial conditions 
leading to non-periodic behavior for larger /i-values (and therefore less potentially interesting 
collective effects). This case essentially shows no dynamical regions structure, as opposed to a 
vast spectrum of dynamical behaviors displayed hy e = 0.9 case. These results suggest that the 
e-values smaller than 0.85 and bigger than 1.0 are far less interesting (at least in the range of 
coupling strength investigated here). 

Visibly, the case of e = 0.9 indeed presents the most interesting example of 2D chaotic stan- 
dard maps network of CCM given by Eq. (j2.3p . both as self-organized complex system and as 
interesting non-symplectic dynamical system. 
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Figure 3.22: Fraction of non-periodic orbits in function of /i-value, averaged over initial conditions. 
CCM on tree is considered with various e-values, as indicated in the picture. 



The CCM with Symplectic CoupHng. The couphng form leading to Eq. ()2.3p is non- 
symplectic by construction, as we favored biological motivation for designing our system over 
preserving symplectic (area-preserving) character of our model for the isolated unit, namely the 
chaotic standard map. Typically, systems of CCM are constructed with non-symplectic coupling 
even for the cases of symplectic units/maps, like logistic map or standard map. However, it is 
to be noted that standard maps can be coupled with preserving their symplectic character, and 
yield a symplectic CCM system of arbitrary dimensionality. Examples include already mentioned 
case studied in [6]: 

and the case of metastable states in weakly chaotic coupled standard maps studied in |12j : 
f x\i]t+i \ _ f x\i][ \ , f Ejec,(yb1t + ^sin(27r2;[i]t)) \ 

and many other. However, these symplectic version of coupling among standard maps generally 
refer only to the cases of clique and chain network structures, without being applicable to general 
network topologies (e.g. scalefree tree). Also, these CCM systems do not exhibit regular behavior 
in the sense of our study here, the phase space dynamics remains chaotic as in Fig. 12. lb (although 
they exhibit other collective effects related to measure-preserving nature of CCM system, like the 
anomalous diffusion [6] or metastable states |12j). 

In order to study the collective dynamics with possible implications for the real complex 
systems, one needs to go beyond the restrictive symplectic coupling of 2D units. For these reasons 
we have not devoted a special attention to the symplectic system of standard maps, and we shall 
continue our examination for the context of previously defined non-symplectic system of CCM. 



Chapter 4 

Dynamic Stability of CCM on 
Networks 

We study the stability of the emergent motion of our network of 
CCM using three different techniques. We characterize the dynamical 
regions further, by observing sharp distinctions between them in 
terms of stability patterns. A variety of dynamical phenomena related 
to non-symplectic coupling are investigated, including different types 
of strange (nonchaotic) attractors and quasi-periodic orbits. Edges of 
dynamical regions are explained through parametric instability. 

We investigate the stability of the emergent dynamics of CCM Eq. ()2.3p on non-directed 4- 
star and tree, considering standard approaches used in nonlinear dynamics |119j . We are dealing 
with chaotic maps that due to mutual interactions achieve regular behavior in a wide range 
of coupling strengths, as shown in detail in the previous Chapter. However, it is of interest 
to see whether this regular (as opposed to chaotic) motion is also stable with respect to small 
perturbations. By characterizing the stability of our system, we establish a further distinction 
with the case of the uncoupled maps which is known to be strongly chaotic. Quantitative results 
indicating emergent motion of nodes to be stable in the sense of nonlinear dynamics, would prove 
the ability of the considered networks to stabilize the dynamics of very chaotic isolated units. 
It is important to note that biological complex systems are generally stable and robust to the 
environmental perturbations, which is necessary for their successful functioning. This implies 
any complex dynamical system modeling a real biological system ought to exhibit a stable and 
coherent behavior. 

We will examine the stability of our CCM system mainly focusing on 4-star which, as demon- 
strated already, contains the core of scalefree tree's dynamical properties. In particular, we will 
employ the following methods: 

1. Finite-time Maximal Lyapunov Exponents (FTMLE), used for general characterization of 
dynamics for all the network structures, based on direct computation of divergence of nearby 
trajectories for a single node 

2. Standard Maximal Lyapunov Exponents (SMLE), computed by the usual technique of Jaco- 
bian matrix time-evolution, used for precise measurement of Maximal Lyapunov Exponents 
for 4-star considered as 8-dimensional dynamical system 
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3. The Parametric Instability, explored using eigenvalues of the system's Jacobian matrix at the 
thresholds separating different types of motion, used to investigate generation of instabilities 
through appearance of various types of bifurcations 

By studying the stability of our system's collective dynamics we shall further characterize the 
dynamical regions already defined in previous Chapter, and explain the differences between them. 
Study of stability is related to concrete dynamical phenomena exhibited by the system at a single 
final state after the transient, in contrast to the global statistical properties examined so far. 
The investigations to follow will complement our previous results and complete the picture of our 
system's collective dynamics. 

4.1 Finite-time Maximal Lyapunov Exponents 

Lyapunov Exponent is generally understood as a measure of divergence of nearby trajectories 
[1191 166^ I38j. computed by systematically recording the distance between two initially very close 
trajectories as they evolve. For a specific phase space point and its trajectory related to a given 
dynamical system, the Maximal Lyapunov Exponent (MLE) is defined to be the fastest divergence 
rate from this trajectory found in its neighborhood. MLE tells how fast are nearby trajectories 
diverging from each other or converging towards each other. Divergence of nearby trajectories 
(which translates into positivity of MLE) is usually related to chaotic dynamics, as opposed to 
regular dynamics which exhibits only non-positive MLE (implying nearby trajectories are not 
diverging from each other). 

MLE are computed by different techniques depending on the system's properties: in this 
Section we will focus on Finite-time Maximal Lyapunov Exponents (FTMLE) that are defined as 
the initial divergence rate of the nearby trajectories [33]. We are seeking a method of quantifying 
trajectory divergence that can be applied for any of the 4-star's or tree's nodes, as this will allow 
a comparative study of the stability of these structures. We will consider the usual 4-star's branch 
node in relation to an outer tree node as done previously. For this purpose we examine in Fig. 14. II 
the divergence between three pairs of nearby trajectories exhibited by 4-star's branch node, picked 
from three different dynamical regions (similar curves are found for the case of tree nodes). Each 
curve reports the time-evolution of the distance between two trajectories dt, divided by their 
initial distance do- For the periodic orbit, the distance ratio shrinks with time, while for the 
chaotic and the attractor orbit it initially diverges, but eventually settles around a certain value. 

For Finite-time Maximal Lyapunov Exponent related to some phase space point belonging to 
the node [i] we shall understand the initial slope of the divergence curve, as pictured in Fig. l4.1l 
and following [61]. More precisely, FTMLE A^^^(xo) associated with the point xq = {xo,yQ)[i] is 
defined by: 



where J\f stays for a small neighborhood around the point xq, d{., .) denotes the distance, and Ut 
stands for the discrete time-evolution dynamics given by the system of CCM. Note however that 
the length of the slope should be appropriately selected for each trajectory separately, depending 
on the dynamical region in question. We shall use FTMLE to characterize the stability of motion 
of selected nodes separately, comparing 4-star's branch to tree's outer node. FTMLE represent 
a stability measure related to the initial divergence of trajectories, and are complementary to 
the Standard Maximal Lyapunov Exponents (SMLE), examined in the next Section. FTMLE 
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Figure 4.1: Time-evolution of distance ratio lii(^) for three characteristic points in three dynam- 
ical regions of the 4-star's branch node. 

depend on the node, as opposed to SMLE which consider the entire system. They are suitable for 
the present study as they allow a parallel examination of both structures, referring to the orbits 
exhibited by similarly located nodes (e.g. branch node vs. outer node). 

A single A^^^(xo) does not fully describe the dynamics related to the point xq as it also 
depends on the test-point x € Af. We therefore average over many x € in order to obtain a 
better approximation, according to the following algorithm: 

1. consider a final steady state of a given network, and focus on the state (xo[^],yoH) for the 
chosen node [i] 

2. consider a random point in this point's close neighborhood, with the distance d be- 
tween them given as " Manhattan distance" (for computational simplicity): (i((x, y), (xq, yo)) = 
\x - xo\ + |y - yol < 1 

3. iterate the dynamics for both points systematically recording the rate that their distance 
ratio in Eq. (j4.ip changes in time. Determine the time interval with the clearest slope 
(slope-length) and compute the corresponding slope over this interval (cf. Fig. l4.ip 

4. Repeat the same procedure for few other random points in the same neighborhood, and 
determine the maximal among the obtained slopes. This is our approximation of A^jj^.(xo) 

Furthermore, the exponent A^^2,(xo) computed this way corresponds to only one point xq on the 
examined orbit. Starting with different points on the same orbit, a distribution of A^^^ can be 
constructed to characterize the entire orbit. The average value of this spectrum will be called 
A^Q^ and understood as the characteristic FTMLE for the whole orbit in question: 

■^max ~< ^maxi^o) >xoeorbit ('^•^) 

We will refer to this definition when considering A^q^(xo) for different nodes/initial conditions, 
usually obtained by averaging over many initial conditions. The distribution of A^jj^(xo) and the 
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value of S^^^ a good estimate of the initial stability of the orbit in question. For simplicity, 
we will determine the optimal slope-length for each dynamical region separately (the slope-lengths 
seem to mostly depend on the dynamical regions), and refer to them for all FTMLE computations. 

In Fig. 14. 21 we report the results of evaluation of A^^^(xo) for 4-star's branch node and tree's 
outer node. For each //-value we compute the )4naxi'^o) many initial conditions, and construct 
a 2D color histogram as done previously in Fig. 13.101 (coloration is in log-scale). There is a clear 
equivalence between profiles for 4-star's branch node and tree's outer node, in terms of similar 
distribution patterns (the tree node picture includes less initial conditions as it is numerically more 
demanding). Small /i- values are characterized by large FTMLE 0(1) similarly to the uncoupled 
standard map, whereas the periodic region shows mainly negative FTMLE. In contrast to these 
two dynamical regions, the attractor region exhibits a larger variety of A^^^j-values which can be 
negative and weakly positive O(10~^). These results indicate a clear analogy between the behavior 
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Figure 4.2: 2D color histograms of Xlnax values in function of the coupling fi, obtained considering 
many initial conditions. 4-star's branch node in (a), and a tree's outer node in (b). 



of 4-star's branch node and tree's outer node in terms of dynamical stability patterns. Given the 
network location analogy of the two considered nodes, the result also suggests an equivalence 
between the general stability of the 4-star and the tree. 

As the unstable motion is related to the positivity of FTMLE, it is illustrative to examine the 
fraction of initial condition leading to the orbits with positive FTMLE in function of the coupling 
strength for various nodes/topologies. We construct such plots starting from the data reported in 
Fig. 14. 21 and present the results in Fig. 14. 31 showing the fractions of unstable orbits characterized 
by positive A^^^j-values in function of /x, for 4-star's branch node and the tree's outer node. The 
profiles show similar patterns, specifically for smaller coupling strengths where the overlap is very 
precise. Attractor region shows a qualitative overlap in terms of coupling strength interval with 
a fraction of positive A^^^-values. Quasi-periodic regions in both profiles agree as intervals of 
coupling strength with some positive FTMLE (quasi-periodic orbits with be studied in detail in 
the following Sections). 
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Figure 4.3: Fraction of orbits leading to positive Xl^ax function of the coupling fx for 4-star's 
branch node and a tree's outer node, considered over many initial conditions. 



FTMLE and Dynamical Regions. The fraction of initial conditions leading to positive 
FTMLE describes the properties of motion of the specific node in function of the coupling strength, 
given that the negative FTMLE are basically equivalent to the periodic orbits. We examine the 
relationship between this fraction and the fraction of non-periodic orbits investigated previously 
(Figs. [3171 &: l3.8p . and used to define the dynamical regions. Such comparison establishes a quan- 
titative result concerning the relationship between the dynamical regions and the stability of 
emergent motion within them. The profiles in Fig. 14. 41 compare the results from Fig. 14. 31 and the 
results on non-periodic orbits for cases of 4-star's branch node and tree's outer node. Clearly, 
in case of the 4-star the overlap is almost perfect: essentially all the non-periodic orbits display 
positive FTMLE, while the periodic ones lead to negative FTMLE, as expected. In the case of 
tree the overlap is not precise, but the structure of dynamical regions is evident. 

It appears the stability of motion of tree's outer node is inherently related to the stabil- 
ity of motion of 4-star's branch node, irrespectively of the coupling strength. Moreover, the 
characterization of motion of our network of CCM in terms of dynamical regions seems to be 
confirmed by the stability investigation. Dynamical regions can be equivalently characterized 
through FTMLE: periodic orbits exhibit negative FTMLE, while non-periodic orbits generally 
exhibit positive FTMLE. To be emphasized however, is that positive FTMLE in the region of 
attractors have far smaller values than the ones related to the chaotic orbits for very small fx- 
values (cf. Fig. l4.2p . As expected, this indicates that although unstable and non-periodic, the 
type of motion exhibited in the attractor region by 4-star's branch node is intrinsically different 
from the usual chaotic one. Furthermore, within the attractor region in Fig. 14. 4b (specifically, for 
/i = 0.048) there appears to be a small fraction of orbits that are non-periodic, but display neg- 
ative (or non-positive) FTMLE. This suggest the presence of new dynamical phenomena known 
as strange nonchaotic attractors, and will be examined in detail in the Sections to follow. 
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Figure 4.4: Comparison of the fraction of orbits having positive Xl^ax the fraction of non- 
periodic orbits, in relation to the couphng strength fi, considered over many initial conditions. 
4-star's branch-node in (a), and an outer tree's node in (b). 



4.2 Standard Maximal Lyapunov Exponents 

In this Section we analyze the stability of our system through Standard Maximal Lyapunov 
Exponents (SMLE) computed using the usual technique involving iterations of linearized system 
(Jacobian matrix) [119^ 166). Due to computational restrictions, we will examine only the dynamics 
of 4-star, comparing the results to the ones from the previous Section. As opposed to FTMLE 
computation model involving only a single node, SMLE is associated with the whole 4-star, 
considering it as an 8-dimensional dynamical system. SMLE give a precise characterization of 
coupled system's stability for a given initial condition, whose precision improves with number of 
iterations considered (in contrast to FTMLE where we had to determine to the optimal slope- 
length). An 8D dynamical system has 8 independent Lyapunov Exponents; for the purposes of 
this study we shall limit the discussion to the largest one among them only, called C- 

The following computational algorithm will be employed in order to evaluate the SMLE 
associated with 4-star for a given coupling strength and initial conditions: 

1. select an initial point xq in eight-dimensional phase space of the 4-star, and compute the 
4-star's Jacobian at that point (defining linearization of the system in the neighborhood of 
xo): 

9f(x) 



J(xo) 

with f denoting the full 4-star's map: 



(4.3) 

x=xo 



xt+i = f(xO, x = (x[l,...,4],y[l,...,4]). (4.4) 

2. chose a random 8D unitary vector uq and compute ui = J(xo) • uq. The value 

Ci = In juil 
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represents the lowest order approximation to the SMLE and is recorded. Redefine ui by 
normaUzing it: ui — > 

compute the next iteration of 4-star's dynamics xi = f(xo), the corresponding J(xi), and 
U2 = J(xi) • ui. The value 

C2 = In |u2| 

is again recoded, and the vector U2 is normalized U2 



U2 
|U2| ■ 



4. steps 3 is repeated for the desired number of iterations, improving the approximation of 
the final SMLE. The value 



k=tt 



c 



k=l 



is defined to be the SMLE for the point xq approximated to the t^-order of the 4-star's 
dynamical system given by f . 

The value C = Ci^o) is related solely to the initial conditions xq and independent from the choice 
of uq, at the limit of — > 00. At this limit the value of C becomes the usual Infinite-time MLE 
characterizing a nonlinear dynamical system |119j . 

SMLE will be hence used to fully characterize the dynamics of 4-star in relation to the coupling 
strength. In Fig. 14. 5b , we show the 2D color histogram of SMLE values in function of /U obtained 
over many initial conditions. The profile is complementary to the Fig. 14. 2k providing additional 
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Figure 4.5: 2D color histogram of SMLE ( in function of the coupling over many initial 
conditions for all nodes of the 4-star in (a), and a zoom at the attractor region of positive ^ for 
the profile in (a), shown in (b). 



stability analysis. The structure of three dynamical regions is visible, with values of SMLE 
being similar to the corresponding FTMLE-values. Also, the similarity between named profiles 
confirms FTMLE to be a valid approximation to the real infinite-time SMLE, in the context of 
single-node dynamics. We also report the zoom to the attractor region in Fig. 14. 5b . showing a 
vast spectrum of dynamical behaviors occurring within this coupling interval. For each /x-value 
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there are certain possible motions (depending on the initial conditions), which does not seem to 
be changing continuously with ^u. This indicates 4-star to be extremely sensitive to both coupling 
strength and the initial conditions, exhibiting rich dynamical behavior in relation to both. In 
the next Section, we shall examine the motion for a specific coupling strength value. However, 
in order to reveal the full spectrum of motions, a much finer discretization of /x-values would be 
required for the attractor region. 

In Fig. 14. 61 we report the comparison of non-periodic orbits and positive SMLE for the 4-star, 
averaged over many initial conditions. The profile is almost identical to the Fig. l4.4h . re-confirming 
the general equivalence between non-periodic and unstable orbits for the 4-star. This further 
confirms the relationship between X^^ax C in terms of their correspondence with unstable 
dynamics. The dynamical regions as intervals of coupling strength with distinctive motions are 




1^ 



Figure 4.6: Fraction of orbits leading to positive Standard Maximal Lyapunov Exponents ( for 
all nodes of the 4-star vs. fraction of non-periodic orbits for 4-star, averaged over many initial 
conditions. 

also recognized by SMLE. The positive values of ^-exponents in the attractor region are far smaller 
than the ones related to small ^-values (similarly to the case of FTMLE). The study of 4-star via 
SMLE gave a precise characterization of its dynamics for the considered coupling range, which 
will be used further to investigate specific dynamical phenomena. 

4.3 Emergence of Weakly Chaotic Collective Dynamics 

In this Section we present a detailed analysis of the attractor dynamical region, examining dynam- 
ical patterns exhibited by the 4-star which display weak chaos. We will investigate the properties 
of strange attractors (already introduced in Fig. l3.15t fc d) appearing in this region, from the non- 
linear dynamics viewpoint. For simplicity, we shall focus of a fixed coupling strength of ^ = 0.048, 
which is rather typical for attractor dynamical region. We will also quantitatively characterize 
the quasi-periodic orbits present throughout all dynamical regions. 
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Dynamical Patterns of 4-star with /x = 0.048. In Fig. 14. 7b we show the distribution of 
SMLE C over many initial conditions for the 4-star with fixed coupling strength of = 0.048. 
Three main types of motion can be immediately recognized: periodic orbits having negative C, 
strange attractors displaying positive ^ > 0.1, and another group of orbits with ^ = 0.038. In 
Fig. 14. 7b we show a zoom to the C ^ 0.1 part of the distribution clarifying between different 
attractor types having similar (^-exponents. Each peak in the distribution of SMLE C represents a 





0.1 0.11 0.12 0.13 0.14 0.15 

(b) 

Figure 4.7: Distribution of SMLE after 10^ iterations for many initial conditions for all nodes 
on the 4-star with fj, = 0.048 in (a), and the zoom to the C ^ 0.1 part of this distribution in (b). 



group of initial conditions which lead to a certain dynamical pattern on the 4-star. Other /u- values 
in the attractor region show similar profiles. 

An attractor corresponding to one of the peaks in Fig. 14. 7b (C — 0.144) is shown in Fig. l4.8t the 
orbits exhibited by the hub and two branch nodes are shown, while the orbit of the third branch 
node overlaps with one of the first two branch nodes (depending on the initial conditions). Again, 
as the orbits have two parts, we are showing only a half of each orbit for clarity (see Fig. 13. lb for 
illustration). Other peaks in Fig. 14. 7b correspond to the similar situations that include all three 
branch nodes having attractors as either Fig. 14. 8b or Fig. 14. 8b . while the orbit of the hub node 
always resembles the one shown in Fig. 14. 8b . The strange attractor in Fig. 14. 8b is a proper fractal 
with fractal dimension df = 1.4, whereas orbits of other nodes (specially the hub) do not display 
this kind of phase space organization. This is a clear example of strange chaotic attractor as it 
possesses the usual properties that define it as such |119| . 

The dynamical situations having their (^-exponents in the range shown in Fig. 14. 7b are in 
opposition with the dynamical pattern occurring for = 0.038 (central peak in Fig. 14. 7b ). Here, 
all four nodes exhibit strange attractors (with similar fractal dimensions of df = 1.4) that are 
typically structured as shown in Fig. 14. 91 Depending on initial conditions leading to this dynamical 
pattern, two branch node's orbits overlap as either Fig. 14. 9b or Fig. 14. 9b : regardless of this, the 
C-exponent always converges to C — 0.038 as — > oo. This situation displays much more 
dynamical organization in terms of inter-node correlation of motion, and it is endemic to the 4- 
star (as opposed to previously examined strange attractor in Fig. 14. 8b which appears on the outer 
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Figure 4.8: Attractors exhibited by 4-star's nodes with ( = 0.144. Hub node in (a), and two 
branch nodes in (b) and (c). The fourth node's orbit is identical to either (b) or (c), depending 
on the initial conditions. 
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Figure 4.9: Attractors exhibited by 4-star's nodes with ^ = 0.038. Hub node in (a), and two 
branch nodes in (b) and (c). The fourth node's orbits is identical to either (b) or (c), depending 
on the initial conditions. 



nodes of the tree as well). Although all nodes are exhibiting fractal attractors, the corresponding 
SMLE is still relatively small (cf. Fig. l4.7h ). Similar occurrences of all nodes displaying strange 
attractors are found throughout attractor dynamical region for other /i-values as well, with C,- 
exponent always remaining C, < 0.04. 

These two dynamical patterns were already investigated in terms of their return times dis- 
tributions, confirming them to be examples of dynamical self-organization (see Fig. l3.15T 3). Here 
we add a convergence study of the SMLE related to the dynamical patterns from Figs. l4.8l &: l4.91 
In Fig. 14.101 we show the distribution of (^-exponents in function of the evaluation time t^. As 
expected, in the case of attractor Fig. 14. 81 the convergence pattern is quite irregular with final 
C-value clearing only after ~ O(IO^) iterations, whereas for the attractor Fig. 14. 91 the conver- 
gence is much more smooth and fast. Although unstable, the dynamical pattern described in 
Fig. l4.9l is far more dynamically structured, and merits further investigation as a clear example of 
self-organization in 4-star. 



4.3. EMERGENCE OF WEAKLY CHAOTIC COLLECTIVE DYNAMICS 



61 



0.35 
0.3 
0.25 



0.15 
0.1 
0.05 



435 



10000 
30000 
100000 
300000 



0.144 



0.1445 



0.145 



0.1455 



(a) 



CL 




(b) 



0.02 0.04 0.06 0.08 



Figure 4.10: Distribution of Standard Maximal Lyapunov Exponents ^ for all 4-star's nodes and 
different evaluation times t^, starting from various points on the attractor. Dynamical situation 
from Fig. l4.8l in (a), and from Fig. l4.9l in (b). 



Additional Analysis of Dynamical Pattern from Fig. 14.91 In addition to return times 
study in Fig. l3.15b . we consider again the single node strange attractors shown in Figs l4.8b &: l4.9b 
that are exhibited by a 4-star's branch node at the coupling strength of /i = 0.048. The corre- 
sponding SMLE are computed above and amount to ^ = 0.144 and C = 0.038, respectively. In 
Fig. 14. Ill we report the distribution of FTMLE A^^^, for the mentioned attractors. The precise 
values (for infinite computation time t^) of SMLE lie within the distribution ranges of FTMLE 
and the distribution for Fig. 14. 9b is much more narrow (similarly to Fig. 14. 10]) . The centers of 
distributions from Fig. l4.11l are respectively equal to Xlnax = 0.118 and Xlnax — —0.005 for both 
attractors. Despite FTMLE being less precisely measured (due to need of slope-length optimiza- 
tion) , they still have a precise dynamical meaning: estimate the initial divergence of orbits related 
to a specific node/trajectory, as opposed to SMLE which gives an overall stability approximation 
for the entire system. In particular, while the FTMLE distribution in Fig. l4.11"k confirms the 
unstable character of the strange attractor from Fig l4.8b . the distribution in Fig. 14. Ill ) points to 
a possibility of strange attractor from Fig. 14. 9b exhibiting some form of stable behavior. This is 
to say that stable behavior might exist in a 2D subspace of 4-star's 8D phase space. The MLE 
distribution having both positive and negative values indicates presence of stable attractor points 
together with unstable attractor points. 

The strange attractors having similar kinds of MLE properties are referred to as strange 
nonchaotic attractors (SNA) |90^ I33j. and are characterized by having non-positive MLE while 
still being geometrically strange (e.g. having fractal phase space organization). Their properties 
and routes to their generation have been extensively studied over the last decade [94:\ [9T] , revealing 
them as a peculiar dynamical phenomena. However, their appearance have so far been related 
exclusively to the dynamical systems with external driving. The attractor in Fig. 14. 9b does possess 
fractal geometrical properties and shows a SMLE of ^ = 0.038, but on a short time-scale seems to 
exhibit a mixed unstable/stable nature with respect to the small perturbations. As this indicates a 
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(a) (b) 

Figure 4.11: Distributions of FTMLE A^^^, for the strange attractors from Figs. 14. 8 b fc l4.9b refer- 
ring to branch-nodes of 4-star, considered over many initial conditions in (a) and (b), respectively. 
The distributions mean values Xl^ax =< ^max >x are Xl^ax — O-l^ and X^^ax = —0.01 for (a) and 
(b), respectively. 



possibility of attractor from Fig. l4.9b having SNA properties, we will employ some known methods 
used to test and describe SNA in order to obtain a better characterization, following [61j . 

We firsty observe that this attractor, in accordance with its fractal structure, does not present 
phase space continuity, as we demonstrate in Fig. 14.12^ by overposing 10 successive trajectory 
iterations to the attractor picture. The apparently uncorrelated successive iterations of the map 
still create a highly structured attractor. 4-star is an endogenously driven system: hub node 
and branch node drive each other, with a time-delayed interaction, as pictured in Fig. 12. 31 The 
attractor from Fig. 14. 9b appears on the 4-star's branch node which is driven by the hub node 
that exhibits attractor shown in Fig. 14. 9b . In Fig. 14.121 ) we report the mutual phase-dependence 
between the two attractors, showing its structure that clearly indicates the absence of a smooth 
curve, which is a typical feature of SNA [86j . Following SNA characterization methods exposed in 
[86], we test the non-differentiability of the attractor from Fig. 14. 9b by considering the properties 
of the time series = dxt[branch]/dxt[hub]. In particular, one can define the sum |33j : 

k=m 

Sm=Yl ^h~^ X exp{A^,,(m - k)} (4.5) 
fc=i 

which examines the named time series taking into account the relationship between positive and 
negative tails in the MLE distribution, as the one shown in Fig. 14.111 3. We numerically compute 
the F^ and show the partial sum series Sm in Fig. l4.12l ;. along with the maximal Sm values 
reported in the inset. Both profiles seem to indicate the partial sum series S^yi grows unbounded 
with the number of terms m, which points to a non-smooth non-differentiable attractor, not 
excluding SNA. Finally, an additional feature that attractor from Fig. 14. 9b shares with the known 
SNA in non-periodically driven maps [86] is the saturation of the maximal distance between the 
attractor points, independently of the time gap between them. In Fig. l4.12B we show the maximal 
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Figure 4.12: Properties of the attractor from Fig. 14.9b : a scatter of 10 successive points in (a), 
phase of the branch node xt[branch] plotted against the phase of coupled hub node xt[hub] in (b), 
partial sums \Sm\ and their maxima max\Sm\ (inset) in function of the number of terms m in (c), 
the maximal distance on the attractor against the distance between the iterations A. 



distance between two points on the attractor in function of the time gap A separating them in 
terms of iterations. For each time gap A we are looking for the maximal distance that two points 
separated by A can attain. Besides depending on the part of the attractor (recall that Fig. 14. 9b 
is only a part of the total attractor of this node), the profiles are constant, implying that points 
with arbitrary time gap can be found at the opposite sides of the attractor, in agreement with 
SNA test in [86]. 

The SNA are so far exclusively found in externally driven dynamical systems (maps), and are 
understood as a consequence of non-periodic driving \33\ \9U[ I81j . Our 4-star system is however 
not externally driven, but it includes endogenous driving in the form of time-delayed interaction 
between the branch-node and the hub (i.e. between the branch node and the rest of the 4-star). 
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Nevertheless, as shown above, the dynamical pattern from Fig. 14.91 and specifically the strange 
attractor shown in Fig. 14. 9b . display clear SNA properties in accordance with the usual tests 
[331 [86]. Ill particular, these properties are: 

• non-integer fractal dimension df = 1.4 

• positive and negative tails in the distribution of FTMLE A^^^ shown in Fig. 14. lib 

• non-differentiability demonstrated through profile of partial sum series Sm from Eq. (j4.5p 
shown in Fig. l4.12t 

• absence of phase space continuity of the attractor shown in Fig. 14. 12^ 

• saturation of the maximal distance of the points belonging to the attractor shown in 
Fig.Eli 

Our study of SMLE which considers 4-star as 8D dynamical system, showed small positive 
values for this dynamical pattern, which suggests the global dynamics of the 4-star to be weakly 
chaotic. However, all the evidence presented (Fig. 14. 121) for the case of particular branch-node 
exhibiting strange attractor shown in Fig. 14. 9b . indicates that this might represent an example of 
SNA in an endogenously driven extended dynamical system. 

Quasi-periodic Orbits. As suggested in the definition of dynamical regions in Fig. 13. 71 the 
system of CCM Eq. (l2.3p also exhibits quasi-periodic orbits for various ranges of coupling strength 
(see Figs. l3.2fe i&: b for examples). Quasi-periodic motion is characterized by the presence of two 
or more incommensurable frequencies, implying it is not periodic (a single frequency present), 
and not chaotic either (continuous spectrum of frequencies) |119j . A quasi-periodic orbit densely 
fills a closed curve in phase space, as shown in the additional example in Fig. 14. 131 Our system 




(a) (b) (c) 

Figure 4.13: Quasi-periodic orbit exhibited by 4-star at /i = 0.054. Hub node in (a), and two 
branch nodes in (b) and (c). The third branch node's orbit overlaps with orbit of either one of 
branch nodes, depending on the initial conditions. 

exhibits quasi-periodic orbits for both 4-star and tree in given coupling ranges, with each node's 
orbit structured similarly to Figs. l4.13k &:bfcc (we are again showing only a half of the whole 
orbit, which as usually has two parts). 
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A peculiar characteristic of quasi-periodic motions is given by its SMLE ( which converges to 
zero with the computation time t^. We shall use this to study systematically the appearance and 
properties of quasi-periodic motion in our CCM on 4-star. In Fig. l4.14b we report the fraction of 
initial conditions leading to |C| < 10~^ over the full range of coupling strengths, indicating presence 
of quasi-periodic orbits in various coupling regions, predominantly around n 0.027 (cf. Fig. l3.7p . 
Quasi-periodic motion also arises throughout attractor region, although less often. In Fig. l4.14b 
we show the histogram of ^-exponents for a fixed coupling strength n = 0.054. As in the previously 
shown histogram for // = 0.048 (Fig. l4.7b ). there is a clear separation between periodic/stable and 
attractor /unstable orbits, with both of them showing many peaks corresponding to various types 
of orbits. However, in Fig. l4.14t ) there is an additional peak around C = 0) corresponding to the 
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Figure 4.14: Fraction of quasi-periodic orbits exhibited by 4-star (defined as |C| < 10~^), in 
function of coupling strength in (a), and distribution of (^-exponents over many initial conditions 
for 4-star at fixed fi = 0.054 in (b). 



initial conditions leading to the orbits such as the one shown in Fig. 14.131 

In order to prove quasi-periodic character of examined motion, we investigate in Fig. 14.15^ the 
convergence of ^-exponent to zero with the number of SMLE computation time t^. Distributions of 
^-values for different orbit points quickly peak around zero; furthermore, the ^-exponent averaged 
over many orbit points < C > shows power-law convergence to zero in function of computation 
time (shown in inset in Fig. 14. 15^ ). with a slope of 1. This clearly points to C = at the limit 

— > oo. Also, in Fig. l4.15l 3 we show a frequency decomposition (power spectrum) of a quasi- 
periodic orbit, with visible peaks around many w-values confirming multi-frequency character of 
these orbits. 

Similarly to periodic orbits, quasi-periodic orbits are rare in the dynamics on uncoupled stan- 
dard map. Their presence of in the emergent motion of our CCM system is another clear sign 
of cooperative effects in its dynamics, existing at various network scales and coupling strengths. 
As mentioned, the dynamics of CCM on the tree also exhibits quasi-periodic orbits on all nodes 
simultaneously, in the similar range of coupling strengths. 
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Figure 4.15: Distributions of (^-values over quasi-periodic orbit points for various in (a); inset: 
average C- value < C > against t^. The frequency (power) spectrum of a quasi-periodic orbit of a 
single node in (b). 



4.4 The Parametric Instability between Dynamical Regions 

In this Section we investigate the generation of instabilities at the edges between dynamical 
regions, explaining the passage from periodic/stable region to non-periodic/unstable region with 
the change of coupling strength. We consider orbits of the 4-star's branch node, and examine the 
changes they undergo in relation to the small changes of //-value, looking for the bifurcations that 
transform one orbit type into another. This way we study the presence of different types of orbits 
for the same coupling strength corresponding to different fractions of initial conditions, and the 
transformations between their ratios with the change of /x-value. We characterize the bifurcations 
through the eigenvalues of the 4-star's Jacobian matrix Eq. ()4.3p corresponding to the phase space 
point in question. As 4-star is 8D dynamical system Eq. (j4.4p . we obtain 4 complex conjugate 
eigenvalue pairs which determine the type of bifurcation |119j . 

We start with periodic region at = 0.015 and report in Fig. 14. 16^ the evolution of period- 
2 orbit into a period-4 orbit as is slightly changed (each point represents only a half of one 
orbit). The corresponding eigenvalues are computed at the marked point in Fig. 14. 16^ and shown 
in Fig. 14.161 3: the right-most eigenvalue crosses the unitary circle, confirming this to be a period- 
doubling bifurcation. Other nodes of the 4-star simultaneously undergo the same bifurcation. The 
variations of /i-value among the orbits shown in Fig. 14. 16^ are O(10~^) - as already mentioned, our 
system of CCM is extremely sensitive to both initial conditions and /i-values. We suggest this to 
be the general mechanism of generation of orbits of different periodicities in all dynamical regions, 
which is confirmed by the presence of large period values on 4-star and tree (cf. Fig. 13.18]) . More 
complex dynamical patterns might be similarly generated as a fast sequence of period-doubling 
bifurcations. 

In Fig. 14. 17b we investigate the generation of a quasi-periodic orbit at /x = 0.014 developing 
from a period-2 orbit (again, we monitor only a half of orbit, seeing quasi-periodic orbit devel- 
oping from a fixed point). Accordingly, the eigenvalues of the Jacobian matrix at the bifurcation 



4.4. THE PARAMETRIC INSTABILITY BETWEEN DYNAMICAL REGIONS 



67 




(a) (b) 

Figure 4.16: Sequence of orbits of 4-star's branch node (only half is shown) undergoing period 
doubling bifurcation at ;U = 0.015 in (a), and the Jacobian matrix eigenvalues computed at the 
bifurcation point (marked in red) for all 4-star nodes in (b). 



threshold are shown in Fig. 14.171 ): right-most complex pair of eigenvalues crosses the unitary cir- 
cle, indicating this to be a Hopf bifurcation |119] . The variations of value in Fig. 14. 17^ are again 
O(10~^), showing the systems to be equally sensitive here as well. Hopf bifurcations are known 




(a) (b) 

Figure 4.17: Sequence of quasi-periodic orbits developing after Hopf bifurcation of a stable point 
at 4-star's branch node with fj, = 0.014 in (a), the Jacobian matrix eigenvalues computed at the 
bifurcation point (marked in red) in (b). 
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to often exist in dynamical systems of this sort, and they are likely to be the general mechanism 
behind the creation of quasi-periodic orbits at all coupling strengths. 

Finally, we examine the creation of a strange attr actor at // = 0.048 corresponding as the 
one shown in Fig. 14. 9b . arising from a periodic orbit by fast increase of orbit's periodicity. We 




Figure 4.18: Formation of the attractor from Fig. 14. 9b through increase of periodicity of a periodic 
orbit at fx = 0.048. First stage corresponds to = 0.048 — 10^ and the last to = 0.048, while 
other stages are between these two /i-values. 

show the stages of attractor formation in Fig. 14.181 with each stage obtained for the same initial 
conditions, but with ^ti- value that vary from II = 0.048 - 10~^ (first stage) to = 0.048 (last 
stage). The process of attractor formation starts with a periodic orbit and displays a sequence 
of several attractors before it is fully formed. Some intermediate attractors (two middle pictures) 
have fractal dimensions smaller than 1; attractor reaches its proper fractal dimension of = 1.4 
only in the last stage. 

The Jacobian matrix corresponding to the final stage of attractor formation from Fig. 14.181 is 
examined, and the results are shown in Fig. 14.191 The eigenvalues for various attractor points are 
interchangingly displaying a stable character (all eigenvalues inside the unit circle as in Fig. 14. 18^ ) 
and an unstable character (one eigenvalues outside the unit circle as in Fig. 14. 181 3) . without any 
apparent correlation. The same pattern also occurs if the Jacobian matrix's eigenvalues are com- 
puted at any of the previous stages (except for the first one, which is always stable). As opposed 
to the previous cases of bifurcations, in Fig. 14.181 one is unable to identify a precise //-value when 
the attractor is formed; accordingly, the stable/unstable pattern of Jacobian matrix eigenvalues 
persists at smaller /x-values as well. Moreover, the stable combinations of eigenvalues (Fig. 14. 18b ) 
appear far more frequently than the unstable ones, which is even more prominent at the earlier 
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stages of attractor formation, which are closer to the periodic orbit. This as an additional ar- 
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Figure 4.19: Eigenvalues of the Jacobian matrix for some attractor points at the final stage of the 
attractor formation shown in Fig. l4.18l A typical stable situation in (a), and a typical unstable 
situation in (b), which are interchanging from iteration to iteration. Jacobian matrix eigenvalue 
with the largest absolute value represents the first-order approximation of SMLE (t^ = 1). 

gument towards the interesting nature of dynamical phenomena from Fig. 14. 9) the structure of 
eigenvalues reports about the stability of linearized system, and in particular, eigenvalue with 
the largest absolute value corresponds to the one-iteration approximation of SMLE ^ (t^ = 1). 
This indicates the nature of linearized version of our system for the dynamical situation from 
Fig. l4.9l to indeed show both stable and unstable nature, as suggested by the FTMLE distribution 
in Fig.SHJa. 

The examined bifurcations refer to the orbits exhibited by the 4-star's branch node; however, 
the transformations of orbit structure occur simultaneously on all nodes. The analysis presented 
in this Section accounts for formation of different dynamical phenomena by all nodes and at all 
coupling strengths (with possibly other types of bifurcations also playing a role) , hence explaining 
the evolution of dynamical regions that network of CCM undergoes with variation of coupling. 
Due to the already established relationship between tree's and 4-star's dynamics, we suggest 
the same mechanism to be behind the transformations of collective motion between equivalent 
dynamical regions for the tree, as the bifurcations seem to be unrelated to the topology details. 
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Chapter 5 



Dynamics on Directed E.Coli Gene 
Regulatory Network 

Additional study of our system of CCM is done using the directed 
gene regulatory network of Escherichia Coli (E.Coh), introduced 
in Chapter [2j Statistical examinations are performed, revealing 
differences and analogies with non-directed topology. Same network 
is also employed for investigation of two-dimensional Hill model 
of E.Coli's gene regulatory network. We particularly examine the 
flexibility and stability of the emergent behavior, along with its 
robustness to fluctuations of the environmental inputs. 



In this Chapter we will investigate the collective dynamics on the largest connected component 
of the directed gene regulatory network of E.Coli shown in Fig. 12. 51 [1001 [70] . with N = 328 
nodes. Two dynamical models will be used: system of coupled chaotic standard maps Eq. (j2.3p 
examined so far, and the continuous-time 2D Hill model |118j of gene interaction introduced in 
Eqs. (|1.6p &: (jl.7|) . We will examine the dynamics of both models mainly employing the techniques 
used above. We are comparing two two-dimensional models with very different dynamical natures: 
while the system Eq. (l2.3p involves strongly chaotic maps, the system Eqs. (ll.6]l &: (ll.7|) displays only 
regular dynamics regardless of coupling. As opposed to investigations done in previous Chapters, 
E.Coli network is directed, meaning that interactions among the nodes are generally not mutual, 
which will affect the properties of the collective dynamics. This network topology also includes 
topological self-loops, allowing the a node/gene to influence itself. It is also to be noted that 
E.Coli network is an example of a real biological network which is empirically found rather than 
computer generated. This implies each node represents a specific gene of E.Coli's gene regulatory 
network, and therefore has a given biological purpose which is encoded in its network location. 

5.1 CCM on Directed Network of E.Coli 

In this Section we consider the dynamical model of CCM given by Eq. ()2.3p on the directed 
network of E.Coli. We adjust the CCM model to the directed nature of the underlying network 
by excluding the degree-normalization factor in coupling form (as some nodes on directed network 
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may not have incoming links, their in-degree is zero). The modified model reads: 

and contains all the key elements of the CCM Eq. (j2.3p (phase coupling, time delay etc.). The prime 
(') refers to the update of the uncoupled standard map Eq. (|2.ip . and the network neighborhood 
ICi includes the nodes [j] = 1,...,N having directed links towards the node [i] (which may 
include the node [i] itself). We will investigate the dynamics of Eq. (j5.ip using the same numerical 
implementation procedure and study approaches from the previous Chapters. Non-periodic orbits, 
time-averaged orbits, node by node FTMLE X^^ax return times will be examined. The studied 
coupling range will include [0, 0.08] as before. Initial conditions are selected randomly from 
{x,y) € [0,1] X [—1,1], and the averaging over them as previously includes 1000 random initial 
conditions. We shall often examine color plots showing certain features for each node; in such 
cases the nodes will be enumerated from to 327. 

We start by reporting the fraction of non-periodic orbits in function of coupling strength /i 
averaged over initial conditions in Fig. l5.li The profile structure at small ^-values resembles the 
profile for non-directed tree: after a quick initial transient the orbits of all nodes for all initial 
conditions stabilize into periodic orbits. However, with increase of coupling strength, a certain 
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Figure 5.1: Fraction of non-periodic nodes/orbits after a transient of 10^ iterations in function 
of the coupling strength averaged over many initial conditions. Inset: zoom to the dynamical 
transition region around fi ~ 0.022. 



fraction of the network nodes destabilizes again into non-periodic orbits (inset in Fig. [5T] provides 
detailed view of this dynamical region with /x ~ 0.022). As the /x- value is increased further, 
the network of CCM re-stabilizes in a short region around fi ~ 0.033, after which remains with 
a relatively constant number of non-periodic nodes for the remaining part of studied coupling 
range. The concrete properties of specific node's orbits will be studied in more detail later. 

In Fig. l5.2h we show the 2D color histogram of y[z]-values in function of the coupling strength, 
equivalent to Fig. l3.10l ;. which confirms the dynamical region structure from Fig. l5.1i For small 
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/i-values the system presents fully clustered organization of periodic orbits in analogy with non- 
directed case, while for larger /i-values a destabilization is visible. Note that certain clusters are 
re-appearing in destabilization region. We also examine FTMLE defined as in Eqs. (j4.i pfc ()4.2p 
by computing A^^^-values node by node, following the same procedure explained in Chapter [H 
The 2D color histogram equivalent to Fig. 14. 21 is reported in Fig. 15. 2b : the same dynamical region 
structure is again visible in relation to the coupling strength. However, as opposed to previously 




(a) (b) 

Figure 5.2: 2D color histogram of -values in function of coupling strength /x for a single initial 
condition in (a), 2D color histogram of A^^^j-values in function of fi, averaged over many initial 
conditions in (b). 



considered non-directed case, after destabilization some nodes show Xl^^x > ^ indicating presence 
of strongly chaotic behavior similar to the uncoupled case. Nevertheless, majority of the nodes 
have Xlnax < pointing to periodic orbits, and some nodes display Af„^^-values in a wide range 
between and 1.5, allowing the possibility of various dynamical behaviors of nodes. 

In the reminder of this Section we perform a detailed analysis of cooperative dynamical pat- 
terns at two particular coupling strength values: at the onset of destabilization at fi = 0.022, and 
the irregular region at /x = 0.05 characterized by Xl^ax > 1- 

Onset of Destabilization at fi = 0.022. We consider the directed dynamics of CCM 
Eq. (|5.ip on the E.Coli network Fig. 12. 51 near fi = 0.022. In Fig. 15. 31 we examine the profile of the 
transition considering the fraction of initial conditions leading to non-periodic orbits for each node 
separately. As the coupling strength increases, the instability develops on a specific group of nodes 
rather than over the whole network. Even after the transition, the instability does not spread to 
the rest of the network, but remains confined to the initially destabilized sub-network. This is in 
sharp contradiction with what observed in the case of non-directed dynamics on tree, where the 
dynamical regions are characterized by all the nodes sharing similar dynamical behaviors. 

We investigate the structure of the unstable sub-network by showing a graphical represen- 
tation of the E.Coli network (from Fig. l2.5p with different coloration for the unstable nodes in 
Fig. 15. 4b . The unstable sub-network includes the hub node (having biggest number of in-links 
and out-links) and other surrounding nodes which are integrated with the rest of the network. 
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node [i] 

Figure 5.3: 2D plot of the instability transition near /i = 0.022: fraction of non-periodic orbits 
for each node (enumerated from to 327) as function of coupling strength averaged over many 
initial conditions. 



In Fig. 15. 4b we show the destabilized sub-network and denote each node/gene with its biological 
name. The unstable sub-network has both incoming and outgoing links with the rest of the net- 




(a) (b) 

Figure 5.4: Instability transition at /i = 0.022. E.Coli's largest connected component with desta- 
bilized nodes (marked in red) in (a), the sub- network undergoing destabilization with nodes/genes 
biological names in (b). 

work, suggesting this to be a non-trivial dynamical phenomenon related to this particular CCM 
system and the architecture of the directed E.Coli network. 
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In Fig. 15. 51 we report three examples of final dynamical steady states at /i = 0.022 by plotting 
10000 iterations of each node on the same phase space (as done in Fig. l3.9p . Besides usual cluster- 
organization involving large majority of nodes, the system also exhibits other kinds of attractors. 
Specifically, the attractor visible in the upper-most region of the phase space corresponds to the 
hub node, and its position shifts depending on the initial conditions, although it qualitatively 
remains the same. To investigate this further, we show in Fig. l5. 61 three typical attractors/orbits 



> 




> 




> 




Figure 5.5: Three examples of final dynamical steady states of all nodes of directed CCM Eq. (j5.ip . 
each corresponding to a single initial condition at fi = 0.022. Shown are 10000 iterations of each 
node on the same plot. 



appearing on various network nodes at this coupling strength. In Fig. l5.6h we see the fractal 
structure of the hub node's attractor, indicating this to be a strange attractor. At this coupling 
strength the hub node always settles in a strange attractor of this type regardless of the initial 
conditions, while no other node ever develops the same attractor. Other network nodes either 
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Figure 5.6: Examples of single-node orbits taken from data in Fig. l5.5i The strange attractor of 
node 2 (hub) in (a), and typical orbits of node 214 in (b) and node 46 (c) (as previously, we do 
not show the entire orbit but only its representative part). 




display usual cluster-oscillations identical to the non-directed case (cf. Fig. 13. lb ), or develop 
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weakly chaotic orbits as the one in Figs. l5.6b fc c. also present in the tree's non-directed dynamics 
(but with larger coupling strength). 

We examine the statistical properties of various single-node orbits by computing their return 
time distributions with respect to phase space partition in x-coordinate. The results are shown 
in Fig. l5.7l where we compare the distributions of return times for the orbits from Fig. l5.6[ All the 




At 

Figure 5.7: Distributions of return times to phase space partition (100000 cells) in x-coordinate 
for single-node orbits from Fig. 15. 61 with fj, = 0.022, fitted with q-exponential distributions. 

distributions indicate presence of long-term correlations in their motion, and can be fitted with 
the q-exponential function Eq. (j3.6p with various values of q. 

Furthermore, we examine the time-averaged orbits y[i] in Fig. 15. 8b by showing the histogram 
of y[i]-values for each node separately. Note that hub node (node 2) has y[i]-value out of usual 
cluster y[«]-values (which are almost integer values), whereas other nodes mostly belong to three 
central clusters with y = — 1, y = and y = 1. FTMLE A^^^, defined as in non-directed case 
are studied for this coupling strength node by node. In Fig. 15. 8b we report 2D color histogram 
showing distribution of A^^^j-values averaged over many initial conditions, for each node. Clearly, 
the destabilized nodes exhibit positive FTMLE (cf. Fig. l5.3p . while other nodes mainly show 
Knax ~ —0.02. Also, unstable nodes show X^^ax ^ 0-3 which is much less than the uncoupled 
case with Xl^ax ~ Moreover, each unstable node shows a specific range of Af^^^j-values, 

indicating a presence of self-organized network behavior at this coupling strength displayed as a 
clear departure from a very chaotic nature of the uncoupled standard map. 

Finally, we study the flexibility or emergent motion by considering the distribution of y[i]- 
values node by node for many initial conditions, and computing the standard deviation cr[i]{y) 
for each node. Small c[«](y) indicates node prefers to settle in a specific cluster only, while large 
(T[i]{y) points to a certain flexibility of node's final motion. In Fig. 15.91 we show the E.Coli's net- 
work, indicating the flexibility of nodes. While centrally located nodes with many ingoing and 
outgoing connections are typically having small iT[z](y) (blue), outer nodes are generally more 
flexible (red) with bigger (T[z](y), i.e. can be found in various motion types depending on the 
initial conditions. This suggests that the network exhibits a self-organization in terms of each 
node showing a certain range of possible behaviors that depend on initial conditions. 
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Figure 5.9: Standard deviation o"[z](y) for each node averaged over many initial conditions for 
fi = 0.022. Red: a[i\{y) > 1.26, yellow 1.15 < a[i\{y) < 1.26, green 1.02 < a[i\{y) < 1.15, blue 
a[i](y) < 1.02. 
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E.Coli Network Dynamics at = 0.05. In addition to the instability transition around 
H = 0.022 we study the instabihty region around fi = 0.05 where some nodes remain clustered with 
periodic orbits, while some other nodes develop fully chaotic motion equivalent to the uncoupled 
dynamics. However, except for full chaos, at this coupling strength the network displays various 
types of non-periodic motion which are specifically dependent on the node in question. 

In Fig. 15. 10^ we show the E.Coli network's largest connected component with red color for 
the non-periodic nodes at /i = 0.05: again, the irregular behavior is localized to a sub-network, 
although this sub-network is fully connected to the rest of the network. In Fig. l5.10T 3 we report 
the graphical representation of the destabilized sub-network, with respective biological names for 
the nodes/genes. In addition to other nodes, this unstable sub-network includes the unstable 




(a) (b) 

Figure 5.10: Instability region at /i = 0.05. E.Coli's largest connected component with destabi- 
lized nodes (marked in red) in (a), the sub- network undergoing destabilization with nodes/genes 
biological names in (b). 

sub-network for fi = 0.022 from Fig. 15. 4b . 

Contrary to what was observed in the case of = 0.022, at this coupling strength the hub 
node does not display any specific attractor, but it exhibits non-localized strongly chaotic motion 
regardless of the initial conditions. However, the behavior of other non-periodic nodes shows a 
variety of structural patterns as illustrated in Fig. 15. Ill where we show three typical node orbits. 
Interestingly, the orbit in Fig. 15.11^ resembles the chaotic attractor found in the non-directed 
4-star case from Fig. 14. 8b . This attractor appears to be inherently related to this particular 
coupling form among the standard maps, rather than to the directedness of graph. The other 4- 
star attractor pattern Fig. 14. 91 was not found in the directed case. We consider the distributions of 
return times to phase space partition in x-coordinate for the single node orbits from Fig. l5.1l1 and 
show the results in Fig. 15. 121 As previously, the distributions are fitted with q-exponential forms, 
indicating presence of similar collective effects at this coupling strength. Note that respective 
values of q in both cases are in the same range, although they always depend on the node. 
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Figure 5.11: Examples of single-node orbits of CCM on E.Coli's directed gene regulatory network 
with = 0.05. Node 18 in (a), node 86 in (b) and node 119 in (c) (as usually, we do not show 
the entire orbits but only their representative parts). 
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Figure 5.12: Distributions of return times to phase space partition in x-coordinate (100000 cells) 
for single-node orbits from Fig. 15. Ill with fi = 0.05, fitted with q-exponential distributions. 



The global statistical properties of the directed CCM Eq. (j5.ip at /x = 0.05 are shown in 
Fig. l5.13l where we show 2D color histograms of y[i]-values and FTMLE A^^^j-values. As opposed 
to /z = 0.022 case, the statistics of y[i]-values depends on the node in question, with different 
nodes exhibiting very different y distributions. It seems each node has a certain "role" in the 
network' collective motion, given by its spectrum of possible behaviors in function of the initial 
conditions. The distribution of A^^^j-values is also node dependent; while most of the nodes 
are stable with Xlnax < 0) some nodes display a double nature by showing (depending on initial 
conditions) either positive or negative FTMLE. In particular, the node 18, similarly to 4-star's 
branch node, can display an attractor (cf. Fig. l5.11"b ) as well as a periodic orbit. Also, while 
some unstable nodes show large FTMLE X^ax ~ ^ similarly to the uncoupled case, some other 
nodes show moderate FTMLE values of Xl^ax ~ O-^ relating them to the /i = 0.022 case. Still, the 
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majority of nodes show a stable behavior with negative FTMLE. We also examine the network 
distribution of standard deviations with data from Fig. 15.13^ . describing the flexibility of 

each node. The results are reported in Fig. 15. 141 where, as previously, we attach a color to each 
node expressing its cr[i] (y)-value. In contrast to /i = 0.022 case, many outer less connected nodes 
are now fixed with very small a[i]{y)-values, while inner more connected nodes display flexibility 
of their final dynamical states in relation to the initial conditions. 

This indicates the unstable motion of specific nodes for the coupling strength of = 0.05 can 
exhibit different patterns that however do share common statistical properties. We illustrate this 
further in Fig. 15.151 where we show different orbits corresponding to various initial conditions for 
three network nodes. The dynamics of a given non-periodic node does maintain certain motion 
properties: note the similarities among the orbit structures and locations in all three cases. Further 
examination of the same nodes reveals constant presence of the given type of motion, which is 
qualitatively independent of initial conditions. Some other nodes show this behavior displaying 
different motion patterns, although most of the nodes are still periodic (cf. Fig. 15. 13b ). This kind 
of motion flexibility is characteristic for biological collective dynamics: e.g., activity of a gene is 
expected to be adaptive to the cellular needs while maintaining the key behavioral properties. We 
furthermore examine the distributions of return times with respect to a phase space partition in 
x-coordinate. The shown profiles refer to the node orbits emerging for different initial conditions 
for node 86 (cf. Fig. 15. 151 3) and node 119 (cf. Fig. 15. 151 ;;). The distribution profiles show minor 
variations, while still maintaining common structural properties regardless of initial conditions. 

This suggests the emergent dynamics of given unstable nodes is not fully irregular, but seems 
to posses certain flexibility with respect to node/initial conditions. As observed already, a feature 
of this type might be related to biological origin of the studied network, in terms of node's/gene's 
"adaptability". As we shall describe in the reminder of this Chapter, the flexibility of emer- 
gent dynamics with respect to initial/boundary conditions is the crucial characteristic of the 
continuous-time Hill model of gene interactions. 
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Figure 5.15: Examples of final single-node orbits (different colors) in relation to different initial 
conditions at // = 0.05. Node 56 in (a), node 86 in (b) and node 119 in (c). 



5.2 Hill Dynamics of E.Coli Gene Regulatory Network 

In this Section we study the two-dimensional Hill model of gene regulation dynamics on the largest 
connected component of E.Coli directed network from Fig. 12. 51 The behavior of a single gene 
(node) will be described as a continuous-time 2D dynamical system defined in Eqs. (|1.6p &: (jl.7p . 
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(a) (b) 

Figure 5.16: Distributions of return times with respect to a phase space partition in x-coordinate 
(100000 cells) for orbits corresponding to various initial conditions with fi = 0.05. Node 86 in (a) 
and node 119 in (b). 



Besides the directed nature of network, the interactions among the genes (nodes) are given as 
activation or repression, which we show in a new graphical representation of E.Coli network that 
includes interaction types as colors of links Fig. l5.17l (following |100^ 170]). 

Two versions of the system are considered, differing by the logic gate model used for multiply 
regulated genes: the SUM model and AND model will have all their logic gates described by the 
SUM-gates and AND-gates, respectively. The dynamical state of a single node/gene [i] within a 
network with nodes is given by the pairs of variables ((?[i],p[«]), and evolves with time according 
to: 

^ = G^{Fu{p[j]),...,FN^{p[j])}-l?q[i\ 

(5.2) 

dp[i] p p 

—— = af q[i] - -ffp[t] 

where in the case of SUM-gates the operator G is given as the average value of all the inputs: 

j=N 

E 



k^^ being in-degree of the node [i] , while in the case of AND-gates G is defined to be the minimal 
of all the inputs [5] : 

gf= mm F,Mj])- (5.4) 
j=l,...,N 

The function Fji describes the way node [j] influences node [i] in case a link between the two 



5.2. HILL DYNAMICS OF E.COLI GENE REGULATORY NETWORK 



83 



nodes exists, and can be either activatory or repressory: 



activation 



(5.5) 



R..F..nji 



repression 



where T denotes the interaction threshold, n is the Hih exponent ( |118| ). ^ allows leaky transcrip- 
tion and (3 defines the maximum expression level (cf. Eq. (ll.7p ). 

A SUM-gate corresponds to the usual paradigm of coupled maps systems as it involves aver- 
aging all the transcription inputs, whereas an AND-gate models the fact that a gene's expression 
rate generally depends on the " worst" transcription input (minimal activation or maximal repres- 
sion). As already mentioned, real logic gates involved in E.Coli network are typically complicated 
and depend on the gene [5] , although can sometimes be seen as a combination of SUM and AND 
gates. For the purposes of this work, we will assume all network gates to belong to only one of 
these categories. 

Despite being non- linear, the dynamical system given by Eq. (j5.5p is fully regular for all the 
parameter value options. We are therefore expecting the emergent dynamics of the coupled sys- 
tem Eq. (j5.2p to be regular as well, which is required for description of a biological system. For 
this reason we will include the arguments like network's flexibility and robustness to noise in the 
reminder of this Section. 

The Numerical Set-up of the Model. We consider the largest connected component of 
E.Coli directed interaction network with N = 328 nodes (Fig. l5.l7|) with the state of each gene 
given by {q[i\{t) , p[i]{t)) at time t. The system's time-evolution is governed by the Eq. (l5.2p under 
the network and parameter constraints. In order to model the cellular input to which the gene 
regulatory network responds, we add a regulatory interaction by an external transcription factor 
(ETF) to every gene that does not have any incoming links (cf. Fig. ll.2p . or has only one coming 
from a self- loop (dark green nodes in Fig. 15.17]) . An ETF for a given gene consists of a fixed value 
of p which is influencing the corresponding gene through an interaction of a fixed type. This way 
we are able to monitor the response/adjustment of the network in relation to the environmental 
inputs defined by the sequence of ETFs, regulating those genes that do not have incoming links 
from the rest of the network. 

The equations are integrated using classical 4*''-order Runge-Kutta integration method, in its 
time dependent version. Given the non-stiff nature of our equations, this explicit method provides 
sufficient precision of the results. For simplicity, we fix the values a = 7^ = 7*^ = 1 for all genes, 
and quench the parameters describing the interaction properties by selecting randomly for all the 
links (both gene— >gene and ETF^gene): 

• /?ji-values from log-uniform distribution on the interval [1, 10] 

• Tjj-values from log-uniform distribution on the interval [1, 10] 

• njj-values between {2, 3, 4} (only cooperative binding is considered) 

• .^ji-values from log- uniform distribution on the interval [0.1, 1] 

We employ log-uniform scale (uniform distribution on the log-scale) in order to emphasize the 
small values of the variables, i.e. the probability for bigger values exponentially decreases. The 
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Figure 5.17: Largest connected component of the E.Coli gene regulatory network (cf. Fig. [23]) 
with the interaction types pictured as red (activatory) or blue (repressory). The nodes without 
incoming links (expect for self-loops) are in dark green, and other nodes are in light green. 



initial states of genes and the values of ETFs are the only quantities that vary from simulation 
to simulation. The computation algorithm used to run a simulation is the following: 

1. upload the adjacency matrix describing the network and quench the (pre-selected) values of 
the parameters as described above 

2. set the values of ETFs by choosing them randomly from log-uniform distribution on [1, 10] 
for each gene without other external input 

3. set the initial conditions on all the genes = 0),p[i](t = 0)) by choosing randomly from 
log-uniform distribution on [1, 10] for both q and p coordinates separately 

4. run the network dynamics as described above until the final state at time tf is reached. 
The considered Runge-Kutta time-step is h = 0.01 (selected as optimal), and a simulation 
involves integration steps. Typically, a simulation is run until tf = 100 

5. consider the final network state {giijit f) , p[i]{t f)) (generally corresponding to the system's 
attractor) and perform the desired data analysis 

In the reminder of this Section we will be employing the algorithm described above for gaining 
insights into the dynamics behind the collective behavior of E.Coli gene regulatory network. For 
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simplicity, we represent the results only using the coordinate p (concentration of protein) which 
is biologically more relevant. The dynamics in (^-coordinate is qualitatively similar. 

Note that similarly to the system of CCM, here we study two-dimensional units on each node, 
running the coupled dynamics in 2 x 328-dimensional phase space. Moreover, each node's local 
dynamics is given by two equations, with equation in q influencing the equation in p (dynamical 
self- loop). This is in analogy with the case of CCM, where two-dimensional standard map also 
involved a dynamical self-loop. 

The Homeostasis and Flexibility of Network's Response. In Fig. l5.18l we show three 
trajectories of 20 representative nodes with fixed ETF- values, but starting from three different 
initial condition on the nodes (values of p[i] and q[i]). After a short transient time, all the 
trajectories settle into a set of constant values which are independent from the initial conditions, 
and related solely to the nodes. We examine only the SUM model here, as the trajectories for AND 
model are qualitatively similar. Repeating this simulation for many different initial conditions, 
we concluded that once the ETF-values are fixed, trajectory of each network node eventually 
settles into a given constant values of p regardless of initial conditions. The set of final values of p 
define the system's attractor. In the context of network-averaged orbit used in previous Chapters 

10 I . : . . . . . . 1 

9 - : 

8 - 
7 

^ 6 - 

-I — ' 

^ 4 

3 - 
2 - ■ 
1 

I ' ' ' ' ' ' ' ' 

2 4 6 8 10 12 14 16 18 

t 

Figure 5.18: Three trajectories of 20 representative nodes of E.Coli network for three different 
initial conditions, with the same ETF-values for the SUM model. 

(cf. Eq. (j3.2|) ). the final system's attractor is therefore a single point, whose domain of attraction 
is the entire 2 x 328-dimensional phase space of all the nodes. We note this to be relevant for 
the biological interpretation: regardless of initial concentrations of proteins and mRNA, all the 
gene expression rates quickly reach their values that correspond to the environmental inputs, 
modeled through EFT- values. The system's attractor is determined exclusively by the sequence 
of ETF-values, and represents system's response to it. 

We furthermore consider the response of the system to a sudden change of ETF-values. In 
Fig. l5.19l we consider the evolution of trajectories of 20 representative nodes, with network's ETF- 
values changed after it settles in the respective attractor (at ti = 10), and then changed again 
after the dynamics settles in a new attractor (t2 = 20). The system responds quickly, with all the 



86 CHAPTER 5. DYNAMICS ON DIRECTED E.COLI GENE REGULATORY NETWORK 
values p[i] reaching new attractor after a short transient time. This is another consequence of the 




5 10 15 20 25 30 

t 



Figure 5.19: Trajectories of 20 representative nodes undergoing two changes of EFT-values for 
SUM model. ETFs are changed instantaneously to a new random set of values at ti = 10 and 
then again at ^2 = 20. 

biological motivation behind the problem: change of environmental conditions induces a response 
from the network, which results in a fast adjustment to the new input. Again, the new attractor 
is given as a constant value of p for each node that is determined solely by ETF- values. 

A question of biological importance regards the network's adaptability to environmental 
changes. In the context of our model, that refers to the range of possible attractors that sys- 
tem can display in relation to various environmental inputs, i.e. ETF-values. To this end we 
examine the distributions of attractor values of p-coordinate for all the nodes for many random 
sequences of ETFs. Results are reported in Fig. 15.201 where we show the 2D color histogram of 
final values of p for each network node, considered over many combinations of ETFs, for SUM 
model and AND model separately. There is a considerable difference among various nodes in 
both models: while some nodes are relatively fixed to a specific value of p-coordinate, other nodes 
are more flexible and exhibit a vast range of possible attracting values of p. The SUM model 
is visibly more adaptable than AND model, as it shows far less nodes with inflexible behavior. 
This is expected, as a SUM-gate includes all the node's inputs, while an AND-gate considers the 
smallest one only, which easily leads to final p = situation. For illustration of nodes' responses 
to various ETFs, we show in Fig. 15. 211 the distributions of attracting values p[i] for four typical 
nodes in SUM model and AND model. Distributions mostly have sharp edges and prominent 
peaks, but always show at lest some flexibility to ETF-values, which in some cases extends to a 
vast range of values of p. Two profiles for a same node in two models have generally nothing in 
common, with some exceptions like the node 89 in Fig. 15. 211 

To complete the investigation of network's flexibility, we study the transient times needed 
for all trajectories to settle in their attractors, once a new sequence of ETFs is applied. We 
run the dynamics for a very long time tf, record the final node values p[i]{tf), and consider the 
time-evolution of the distance: 

d[m = \p[m-pmtf)\ (5.6) 

which measures the speed of [i]-node's trajectory approaching its final value p[i]{tf). In order 
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Figure 5.20: 2D color histograms of attractor values of p for SUM model in (a) and AND model 
in (b), obtained by considering attracting values of p[i] for all the nodes in relation to many ETF 
combinations. 
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Figure 5.21: Distributions of attracting values p[i] (same data from Fig. 15. 201) in relation to ETFs 
for four nodes displaying typical profiles. SUM model in (a) and AND model in (b). 



to illustrate the global network's approach to attractor, we consider the evolution averaged over 
all the network's nodes < >• In Fig. 15. 22^ we report an example of < > evolution 

for a SUM model trajectory and an AND model trajectory using the same ETFs and initial 
conditions. Visibly, the AND model's trajectory is faster to arrive to the respective attractor. To 
gain quantitative insight into the distribution of transient times, we re-define the attractor time 
to be t J such that < d[i]{tf) > < 10~^^. In Fig. 15.221 3 the distributions of transient times tf is 
shown for both models, considered over many ETF-values. The AND model clearly shows a much 
faster approach to the attracting values p[i] regardless of the environmental conditions. 
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t tf 

(a) (b) 

Figure 5.22: The time-evolution of node-averaged distance to attractor < d[i](t) > for both models 
under the same ETFs and initial conditions in (a). Distributions of tj-values with < d[i]{tf) > < 
-|^Q-15 Jqj. ]-,q^]2 models taken over many ETFs and initial conditions in (b). 



It is interesting to note that each model has one of the biologically needed advantages: while 
the SUM model provides more flexibility for network adjustment to changes in the environmental 
inputs, the AND model is much faster to adopt to the changes. This suggests the real biological 
logic gates in the E.Coli gene regulatory network might have characteristics of both of the exposed 
models. 

Stability of the Systems' Attractors. In order to quantitatively characterize the stability 
of an attractor corresponding to a sequence of ETF-values, we examine the divergence of the 
nearby trajectories once the system relaxes in the attracting constant values of p. We compute 
the Finite-time Maximal Lyapunov exponents FTMLE A^^^ (as done previously with network 
of CCM Eq. (|2.3p and Eq. (|5.ip ) for each systems orbit, obtaining a measure of systems' attrac- 
tor stability. The procedure is equivalent to the one exposed previously (cf. Chapter [H) and 
involves time-evolution of the distance between the attractor point and a point in its close 

neighborhood p[i]{t): 

M = d{p[i]{t>tf),p[i]it>tj)) 

do\i] d{pmtf),pmtf)) ^ ■ ' 

thus measuring the trajectory divergence for the node [i]. The FTMLE A^^^. is then defined as 
the initial slope of the divergence curve (on the log scale). The values A^^^, are here equivalent to 
previously considered A^„^. The averaging over the orbit is redundant, as the orbit (trajectory) 
is composed of a single point (attracting value of p[i]). 

In Fig. 15.231 we show the evolution of for a few representative nodes for E.Coli's network 
in their attractor states. All the nodes initially exhibit exponential convergence of nearby tra- 
jectories, after which the distance settles to a fixed value which is smaller than O(10~^^) (and 
therefore small enough to be confused with zero in the context of finite numerical precision) . The 
initial slope of the divergence curve is defined as the FTMLE A^„^, which are measured accord- 
ingly for each node separately. As the attractor properties are solely given by the ETF-values 



5.2. HILL DYNAMICS OF E.COLI GENE REGULATORY NETWORK 



89 



5 


o 

~] -10 
O -15 

-20 
-25 

5 10 15 20 25 30 35 40 45 50 

t 

Figure 5.23: Divergence of nearby trajectories at the attracting state for a few representative 
nodes in the E.CoH SUM model (AND model shows a similar behavior). 



(and the SUM/ AND model in question), the respective values of A^^^, node by node are also 
defined by them. We therefore compute the distributions of A^^^-values for each node, consid- 
ered over many combinations of ETFs, for both models. Results are reported in Fig. l5.24l in form 
of 2D color histograms. Clearly, both models exhibit non-positive FTMLE for all nodes and for 
essentially all ETFs. This implies the finial attractor system's state to be stable regardless of 




ETF-values. This is of large importance in the context of biology: for each functional network 
mode (i.e., for all possible combinations of ETF-values), the network's final steady state is stable 
to small perturbations. This means the network's function is robust to small changes of protein 
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concentrations that may occur due to e.g. thermal noise within the cell. Notably, AND model 
shows somewhat more uniform stability pattern which is less variable from node to node, which 
is a further advantage of the AND model. 

Models' Robustness to the Noise. We extend our model by assuming a weak time- 
dependence of ETF- values. We modify ETFs by adding a gaussian noise to each of them centered 
at pre-selected ETF-value and having intensity (standard deviation) rj. The value of r] is equal 
for all the ETFs, and is much smaller than the values of ETFs themselves (typically O(10~^), 
while ETF- values are O(IO^)). In this way we are closer to the model of the living cell, where 
the presence of external transcription factors giving input to the gene regulatory network always 
fluctuates within some standard deviation. 

As as result of model modification, the system never reaches the precise attracting (constant) 
values of p, but the trajectories p[i](t) of all the nodes [i] always fluctuate. However, the time- 
average of each node-trajectory p[i] considered over a long time interval (and after transient time) 
coincides with the non- fluctuating attracting values of p for the same sequence of ETFs. This 
indicates the network to be robust in the response to the fluctuations (noise), which is crucial 
for the functional operation of biological networks. We illustrate this in Fig. 15. 25^ showing the 
trajectories of fluctuating and non-fluctuating system for the same ETFs as they evolve, and in 
Fig. 15. 251 3 where we examine the fluctuation properties of trajectories for the same two systems 
away from the transients. The noisy system (blue in Fig. l5.25ll settles to fluctuate around the 
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(a) 

Figure 5.25: Fluctuating trajectories of few representative nodes with noise rj = 0.2, in comparison 
with the corresponding non-fluctuating trajectories for SUM model. Time-evolution of trajectories 
towards the attracting values in (a), and the fluctuations of trajectories after transients in (b). 

attracting values of the non-fluctuating system with a speed which is in agreement with transient 
times examined in Fig. 15.221 To confirm this quantitatively, we show in Fig. 15.261 the attracting 
non-fluctuating values of p and the means of p[i] of the fluctuating trajectories for each node, 
which show a clear overlap for all the network nodes, as suggested. Biologically, this indicates the 
system to be able to "on average" maintain the protein concentrations at their respective levels 
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Figure 5.26: The attracting values of p of a non-fluctuating system and the means of p (averaged 
over a long time interval away from transients) for the system with the same ETFs fluctuating 
with 1] = 0.2. 



defined by the sequence of ETF-values, despite the ETF fluctuations. Robustness to noise of this 
type is generally present in all biological systems. 

After transients, all the nodes fluctuate around their corresponding non-fluctuating values of 
p, with the properties that depend on the node in question and the noise intensity rj. We quantify 
the trajectory fluctuations after transients by evaluating their standard deviation crp[i]. The value 
of Cpli] (if computed over a sufficiently long time interval) depends on the node, ETFs and the 
noise intensity rj. In order to systematically study the systems' response to ETF fluctuations, we 
consider the distribution of crp[i] for each node over many combinations of ETFs. The results are 
reported in Fig. 15. 271 for SUM model and AND model separately for the noise of rj = 0.01. While 
the value of crp[i] depends on both node in question and the ETF- value, it is immediately clear 
that E.Coli gene regulatory network in both models acts as a noise inhibitor. While the input 
fluctuations of the ETFs are t] = 10~^, in both models the fluctuations on the node trajectories 
are typically smaller by few orders of magnitude. The AND model inhibits the noise better with 
some nodes output fluctuations being only O(10~^^ (which is partially due to some of these nodes 
always displaying zero attracting value of p[i]). Both models exhibit a specific range of crp[i] for 
each node, which can be attributed to the biological network preferring some nodes more robust 
than the others. In computational terms, the inhibition of ETF fluctuations is due to the shape of 
the Hill's functions: the range of fluctuations of ETFs values p is generally larger than the range 
of fluctuations of F+{p[ETF]) and F-{p[ETF]). 

To illustrate more directly the inhibition of fluctuations which characterizes our models, we 
compute < (Tp[i] >, the average of cTp[i] over all the nodes for a few ry- values, and present the 
results in Fig. 15.281 The value of < crp[i] > is constantly smaller by about 2 orders of magnitude, 
as already observed. Moreover, there appears to be a power-law relationship between the values 
of < ap[i] > and the corresponding value of noise rj, as suggested by the power-law fit in Fig. l5.28l 
with the slope of 1.23. These results for the SUM model indicate noise inhibition to be even more 
efficient for the case of AND model. Visibly, the examined models of gene regulation in E.Coli 
posses a high degree of self-organization and are strongly robust to the fluctuations of input ETFs, 
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Figure 5.28: The node-averaged standard deviation of all the trajectories fluctuations < (Tp[i] > 
in function of the noise rj, fitted with slope of 1.23. Consideration is done for the SUM model. 



maintaining their emergent operation only weakly influenced by the presence of the noise. 

The studied models of gene dynamics of E.Coli gene regulatory network correctly predict the 
main features that a biological system of this sort is expected to have: homeostasis, flexibility of 
response to changes in external transcription factors (ETFs), stability of system's final attracting 
state regardless of particular combination of ETFs and the robustness to the fluctuations of ETFs 
and their inhibition. Furthermore, our finding agree with the fact that E.Coli's gene logic gates 
are a complicated composition of SUM-gates and AND-gates [5], by showing that each of these 
models possesses the biologically required characteristics to a certain degree. 



Chapter 6 

Conclusions 



We conclude the present Thesis, summarizing and discussing the 
results reported in the previous Chapters. A list of open questions 
related to the studied topics is also provided. 

We have exposed the research results according to the study directions established in Chapters 
[TJ and [21 We considered networks of 2D chaotic maps coupled with time delay, with a variety 
of directed and non-directed topologies. Our finding indicate that the considered networks are 
able to introduce regularity and stability in the motion of very chaotic maps due to inter-node 
interactions. Moreover, we found a clear dynamical relationship between a large network structure 
and its typical dynamical motif, arising as a consequence of topological relationship between the 
two structures. We recognize the presence of three main types of cooperative dynamics: regular, 
weakly chaotic, and self-organized motion with MLE close to zero. Ways of transition from 
periodic to non-periodic orbits were also investigated. In our 4-star system that can be seen as 
8D nonlinear dynamical system we have found evidence of specific dynamical phenomena termed 
strange nonchaotic attractors. 

6.1 Outline of Main Results 

In this Section we give a more systematic presentation of the main results discussed in this Thesis. 

The Regularization of Collective Dynamics of CCM on Networks. For any non-zero 
value of coupling strength /i the dynamics of chaotic two-dimensional standard maps on both 
directed and non-directed network becomes regular, which is for small values of coupling strength 
realized in form of periodic orbits with different periodicities. Due to the inter-node interactions 
each node develops a periodic motion after some transient time, despite very chaotic nature of 
the isolated standard map. The mechanism behind the regularization process involves the inhi- 
bition of chaotic diffusion characteristic for standard map, which with time-evolution increases 
the correlations among the motion of the nodes. The transient time of the regularization process 
depends on the size of the network structure and the coupling strength. For large scalefree tree 
the process typically starts from outer less connected nodes, and moves along the tree branches 
towards the hub node. The final dynamical steady state is characterized by periodic orbits on all 
the nodes, which are stable in the sense of having negative FTMLE and have period values that 
might vary from node to node. Clearly, the scalefree tree topology is able to introduce regularity 
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into the motion of 2D strongly chaotic maps, for even very small values of the coupling strength. 
At larger coupling strengths, the motion of the nodes is not periodic, but nevertheless displays 
other self-organizational features. 

Dynamical Regions of CCM in Relation to the Coupling Strength. With increase of 
coupling strength, the network of CCM does not simply maintain the structure of periodic orbits, 
but instead develops a variety of new collective behaviors. In particular, depending on the value 
of coupling strength within the investigated coupling range [0, 0.08], we recognized four dynamical 
regions: chaotic region (0 < /i < 0.012), periodic region (0.012 < /i < 0.026), quasi-periodic region 
(0.026 ^ /U ^ 0.04) and attractor region (0.04 ^ ^ 0.065). The regions are characterized by 
specific types of collective dynamics that can occur within the respective coupling ranges for given 
initial conditions [61]. The types of motion are also somewhat interwoven, in the sense that pe- 
riodic orbits occur at all coupling strengths, for the appropriate initial conditions; quasi-periodic 
orbits also appear at some specific coupling values outside the quasi-periodic region. Stability 
investigation via Maximal Lyapunov Exponents (MLE) over the whole coupling range confirms 
the dynamical regions: chaotic orbits are highly unstable with large positive MLE, periodic orbits 
are always stable with negative MLE, quasi-periodic orbits are neutrally stable with zero MLE, 
and strange attractors are mostly weakly unstable with small positive MLE (although evidence 
of strange nonchaotic attractors have been found, with very small MLE). Directed network pos- 
sesses differently organized dynamical regions, although the periodic region is still present with 
the same properties in a similar coupling range. The changes between the dynamical regions 
can be explained through the parametric instability of an emergent orbit. In particular, we have 
found and investigated the cases of period-doubling and Hopf bifurcations occurring at the edges 
of dynamical regions and inducing the changes between the types of motion. 

Characterization of Self-organization Properties. The dynamics of the network of CCM for 
directed and non-directed structures in all dynamical regions (with exception of chaotic region) is 
characterized by a variety of self-organization effects, which were studied through their statistical 
and stability properties. For all structures, the regularization of dynamics is accompanied by the 
clusterization of periodic (and some non-periodic) orbits. Each node's periodic orbit belongs to 
one of the clusters defined by a certain y- value of the orbits' time-averages y[i]. The y- values 
considered over the entire network make a discrete set, with almost precise integer spacing in 
y-coordinate. The distribution of network distances between the nodes belonging to the same 
cluster is different from the distribution of topological distances on the network, with prevailing 
distance of 2 links [601 [9]. For non-directed topologies, periodic region shows power-law tails in 
the period distribution. The attractor region is characterized by q-exponential distributions of 
return times with respect to the phase space partitioning. Anomalous diffusion in phase space of 
similar sort was found for all values of the coupling strength in the motion of the tree. Directed 
topologies show destabilization of dynamics at larger coupling strengths, with motion properties 
of some nodes being identical to the case of the isolated maps, i.e. chaotic. The destabilized 
nodes are mutually connected and construct a sub-network which is integrated with the rest of 
the network. Nevertheless, directed graph exhibit variety of non-periodic weakly unstable orbits, 
displaying strange attractors and q-exponential distributions of return times to phase space par- 
titions for particular nodes and coupling strengths. 
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Dynamical Patterns in CCM on 4-star. Various manifestations of non-symplectic coupling 
were found in the dynamics of CCM on 4-star motif at particular values of the coupling strength. 
This include a wide spectrum of strange attractors in the attractor region having positive MLE, 
but with much smaller values (0(0. 1)) with respect to the isolated map case (0(1)), and quasi- 
periodic orbits. A dynamical situation was revealed with all four 4-star nodes exhibiting strange 
attractors, and with a distribution of FTMLE with negative peak value having tails on both 
negative and positive FTMLE values. SMLE for this case was found to be C = 0.038. However, 
considering the attractor on a 4-star 's branch node (coupled to the rest of the 4-star), we found 
all the standard characteristics of a strange nonchaotic attractors [33] present. This is the first 
example that the evidence of SNA was found in an endogenously driven system (while SNA are 
so far known to arise only in externally driven dynamical systems) |61j . Collective effects in the 
motion of the CCM on 4-star are a further example of self-organization in our system which results 
in entirely different dynamical patterns than for the case of uncoupled standard maps. 

Dynamical Relationship of Scalefree tree vs. 4-star. The expected relationship between 
the collective dynamics on a large-scale (scalefree tree) and a small-scale (the typical dynamical 
motif capturing the tree topology, 4-star) was found. The dynamical patterns observed on the 
large scalefree tree seem to be induced by the dynamics exhibited on the 4-star, which suggests 
the known topological relationship between these two structures to extend to the context of dy- 
namical behavior. The identified dynamical regions for the tree and the 4-star overlap in terms 
of the coupling strength range, with same motion types being present in all of them (although 
sometimes for different fractions of initial conditions) . Main statistical properties of the emergent 
motion, along with the stability profile in function of the coupling strength is also shared by these 
two structures. By comparing the clustering properties of 4-star and 4-clique motif, we found the 
4-star to resemble much better the clustering properties observed for the scalefree tree. However, 
the clustering properties of both 4-star and 4-clique seem to be shared with the modular network, 
grown from the scalefree tree with the addition of cliques. 

Other Coupling Forms. CCM with different coupling forms were also studied for comparison. 
Specifically, the equivalent non-delayed system showed similar general properties, although with 
less structured collective effects (e.g. less organized cluster structure). This confirmed the advan- 
tage of the time delay in the node coupling for the context of two-dimensional chaotic maps. The 
systems with different standard map's chaotic parameters e were also examined, showing very 
little collective effects (larger e) or showing periodic orbits as basically the only collective effect 
(smaller e). This study confirmed our choice of e- value (e = 0.9) as optimal for the purposes of in- 
vestigation done here, which was primarily oriented towards exploration of the self-organizational 
effects arising in collective dynamics of coupled 2D chaotic maps. 

Tw^o-dimensional Hill Model of E.Coli Gene Regulatory Network. The continuous-time 
model of gene regulatory network of E.Coli based on Hill two-dimensional approach to the gene 
interactions was implemented on the same directed network. The system exhibits the expected 
collective behavior in the context of biologically motivated functional network. This includes 
homeostasis (final attractor state is independent from the initial states of the genes/nodes), flexi- 
bility of response to environmental inputs (quick adjustment to new external transcription factors 
with wide range of respective attractors) and stability of the attractor (negative FTMLE for all 
nodes regardless of ETF- values). The version of the system with logic-gates modeled as SUM- 
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gates shows wider flexibility of response to ETFs, while the system with AND-gates is faster in 
responding to ETF changes. Both models are robust to fluctuation in the ETF-values, reducing 
the standard deviation of fluctuations from ETFs to the node's attractor trajectories by almost 
two orders of magnitude. Directed E.Coli network does not only produce a largely stable system 
with very chaotic standard maps, but also possesses the needed functionality in the context of 
biological systems. 

6.2 Future Directions 

The results presented in this Thesis are by no means exhaustive in regard of study of two- 
dimensional dynamical units on complex networks. In this Section we outline few directions that 
can be taken in future studies, in order to extend/complete the arguments presented here. 

• CCM with alternative coupling forms (e.g. other versions of diffusive angle-coupling or 
possibly action-coupling) might be designed appropriately in the context of specific 2D 
dynamical system in use. Different chaotic 2D maps can be employed, with varying levels 
of chaoticity. As the time delay appears to be relevant for oscillator-motivated dynamical 
systems, a different structure of time delays might be of interest (for example, time delays 
that vary from link to link according to a prescribed rule). Additionally, disordered systems 
might be of interest, where each node/map has a different level of chaoticity (that is for 
example a system of coupled standard maps with varying e- values). 

• The CCM considered here appears to be relatively robust to the underlying topology: it 
would be interesting to examine how far the robustness persists, i.e. for how diverse spectrum 
of networks this CCM would still maintain its key collective properties (clustering, periodic 
orbits, anomalous diffusion, etc.). More generally, investigations of 2D CCM on various 
networks are still very scarce. An open question revolves around the issue of modularity: 
how much insight on the behavior of a large structure can be gained from studying the 
behavior of its building blocks (modules), i.e. by following a "bottom- up" approach? 

• In our study of CCM on 4-star motif we revealed a large spectrum of collective dynamical 
effects, which include various strange attractors with small Lyapunov Exponents. Some of 
them have properties that are possessed by a special category of attractors referred to in the 
literature as strange nonchaotic attractors (SNA). These dynamical phenomena are so far 
known to arise only in externally driven dynamical systems (whereas our system of 4-star 
is endogenously driven). Our evidence opens a question regarding possible mechanism that 
can lead to SNA in dynamical systems without external driving, where further investigation 
is needed. 

• We have demonstrated some results regarding the collective dynamics of 2D units on a 
directed topology of a real biological system - the gene regulatory network of bacterium 
E.Coli. A wider investigation of dynamics on directed structures is still to be done, specially 
in the context of 2D CCM. Furthermore, 2D description of genes and their interactions might 
help design better models of gene regulatory networks, allowing a vast range of studies and 
applications. Realistic model of gene interactions could also be implemented on computer 
generated networks, and followed by studies of optimal design of functional gene networks. 
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We append two papers containing parts of this Thesis, in particular: 

• Z. Levnajic and B. Tadic: "Self-organization in Trees and Motifs of Two-Dimensional 
Chaotic Maps with Time Delay", Journal of Statistical Mechanics: Theory and Experi- 
ment, P03003, 2008 ([61]) 

• Z. Levnajic and B. Tadic: "Dynamical Patterns in Scalefree Trees of Coupled 2D Chaotic 
Maps", ICCS 2007, Lecture Notes on Computer Science 4488, p. 633-640, 2007 ([60]) 
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